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ABSTRACT 


The  need  for  direct  turbulence  measurements  in  the  upper 
ocean  arise  from  widespread  requirements  to  correctly 
parameterize  momentum  and  scalar  fluxes  across  the  air/ocean 
interface.  Until  recently  these  observations  were  limited  by 
a  lack  of  instrumentation  capable  of  measuring  the  fine- 
structure  velocity  field  down  to  dissipation  scales. 

The  recently  developed  CDV  package  allows  simultaneous 
sub-centimeter  resolution  measurements  of  temperature, 
conductivity,  pressure,  shear  and  the  3  component  velocity 
field.  As  part  of  the  development  of  this  system  a  Monte-Carlo 
simulation  was  used  to  analyze  the  performance  of  several 
spectral  estimators  of  the  mean  acoustic  Doppler  shifts,  from 
which  the  velocity  components  are  derived.  The  selected 
algorithms  were  implemented  on  a  Digital  Signal  Processor 
allowing  real-time  estimation  of  the  velocity,  shear  and 
scalar  quantities. 

To  verify  the  performance  of  the  CDV  package,  it  was 
deployed  off  Wharf  2  in  Monterey  for  24  hours  while 
simultaneously  recording  meteorological  data.  ’.<eynold's 
stresses,  buoyancy  fluxes  and  fine  scale  stratification  were 
characterized  and  the  surface  gravity  waves  fiald  identified. 
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I. 


INTRODUCTION 


Turbulence  is  a  feature  of  flows  and  is  characterized  by 
its  irregularity  or  randomness,  diffusivity,  large  Reynolds 
numbers  (highly  nonlinear)  and  by  three-dimensional  vorticity 
fluctuations.  Turbulent  flows  are  also  dissipative,  requir’ng 
a  continuous  supply  of  energy,  and  limiting  the  smallest 
scales  to  be  far  larger  than  any  molecular  length  scale.  The 
smallest  scale  occurring  in  turbulent  flows  is  called  the 
Kolmogorov  length  scale  [Ref.l].  These  characteristics  make 
any  deterministic  approach  to  the  description  of  turbulent 
flow  impossible,  and  observations  very  difficult.  On  the  other 
hand,  statistical  approaches  to  the  equations  of  motion  lead 
to  an  unresolved  system,  requiring  extra  assumptions  in  order 
to  obtain  closure  [Ref.z]. 

In  order  to  characterize  turbulence,  we  need  to  perform 
measurements  of  mean  quantities  and  of  fast  varying 
quantities,  requiring  high  spatial  resolution,  fast  response 
sensors  and  stable  observation  platforms,  A  classical  approach 
to  separating  mean  quantities  from  turbulent  fluctuation,  is 
to  assume  a  spectral  gap,  so  that  by  using  a  running  mean  we 
are  able  to  isolate  perturbations  from  the  mean  flow  (Reynolds 
averaging)  [Ref.2]. 

A  more  complex  procedure  is  required  if  one  wants  to 
separate  a  linear  wave  field  from  superimposed  turbulence 
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[Ref. 3  and  Ref.4].  One  method  is  just  to  subtract,  from  the 
fast  varying  variables,  the  components  that  are  coherent  to 
the  wave  field.  The  cross  spectral  variables  and  its  coherency 
are  used  to  identify  the  wave  components  and  the  correspondant 
frequency  ranges.  This  method  will  only  demodulate 
superimposed  turbulence  from  the  wave  field  assuming  that  no 
energy  is  being  transferred  (linear  approach),  and  non-linear 
interaction  between  the  two  fields  are  not  identified. 

Turbulence  velocity  observations  are  frequently  limited  by 
a  lack  of  instrumentation  capable  of  measuring  the  fine- 
structure  temperature,  salinity  and  velocity  field,  down  to 
dissipation  scales.  The  CDV  package  recently  developed  at  the 
Naval  Postgraduate  School  allows  simultaneous  sub-centimeter 
resolution  measurements  of  temperature,  conductivity, 
pressure,  shear  and  the  3  component  velocity  field.  These 
concurrent  high  resolution  ocean  microstructure  measurements 
will  aid  our  understanding  of  a  range  of  phenomena  that 
requires  the  determination  of  momentum  or  scalar  fluxes. 
Unresolved  problems  include  layer  stability  characterizations, 
turbulent  closure  problems,  dissipation  and  energy  cascade 
estimations,  observation  of  entrainment  and  turbulence 
development,  air-sea  interaction,  and  the  superposition  and 
coupling  of  turbulence  to  mean  advection  and/or  wave  fields. 

As  part  of  this  work,  real  time  signal  processing 
programs,  programed  in  C  environment  are  developed,  for  the 
CDV  package.  This  work  is  focused  in  creating  and  implementing 
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algorithms  that  estimate  mean  doppler  shifts,  from  which  the 
velocity  components  are  derived;  on  the  real  time  processing 
and  filtering  of  the  analog  variables,  including  temperature, 
conductivity  and  shear;  and  on  the  continuous  real  time 
graphical  monitoring  of  selected  variables.  Also,  using  a 
Monte-Carlo  simulation,  the  performance  of  several  spectral 
estimators  of  the  mean  acoustic  doppler  shifts  are  analyzed  to 
guide  the  real-time  processing  algorithm  choice. 

The  CDV  package,  after  being  tested  and  calibrated,  was 
deployed  on  the  2  and  3  August  1991,  in  shallow  water  off  a 
Wharf  in  Monterey  to  test  its  ability  to  measure  momentum  and 
scalar  fluxes.  It  was  fixed  to  a  2  meter  high  bottom  mounted 
tripod,  with  the  instrument  facing  the  surface  at  depths  of 
the  order  of  0,5  meters,  while  simultaneously  recording 
meteorological  data.  Small  swell,  weak  winds  and  overcast  sky 
were  observed  throughout  the  observation  period. 

A  procedure  is  developed  that  allows  the  estimation  of 
Reynolds  stresses  and  buoyancy  fluxes  based  on  the  measured 
quantities.  Spectral  methods  are  used  to  characterize  which 
part  of  the  observed  scalar  variance  fields  are  turbulent  v.s. 
advected  by  the  surface  gravity  wave  field. 

This  thesis  is  organized  as  follows.  Chapter  II  summarizes 
the  theorical  background  and  motivations  for  microstructure 
observations.  Chapter  III  describes  a  series  of  Monte-Carlo 
simulations  that  study  the  performance  of  several  spectral 
estimators  for  the  mean  acoustic  doppler  shifts  used  to 
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estimate  velocities.  Chapter  IV  describes  some  recent  methods 
and  equipment  that  led  to  the  development  of  the  CDV  package; 
the  CDV  package;  the  signal  processing  developed  algorithms; 
and  some  validation  tests  and  error  analysis  performed  with 
real  data.  Chapter  V  describes  the  Monterey  Wharf  #2 
experiment  presenting  the  Reynolds  stresses  and  ocean  fluxes 
estimation  procedures.  Finally  chapter  VI  presents  some 
conclusions  and  recommendations  for  future  work. 
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II.  MOTIVATION  FOR  MICROSTRUCTURE  OBSERVATIONS 


Microstructure  measurements  are  required  for  oceanography 
problems  such  as  the  study  of  mixing  parameterizations , 
dissipation  and  energy  cascade,  surface  layer  entrainment, 
air-sea  interaction  and  coupling,  upwelling  systems  and  mixing 
layer  models.  This  chapter  will  introduce  some  theorical 
background  to  ocean  turbulence  emphasizing  the  experimental 
needs . 

A.  REYNOLDS  EQUATIONS 

The  dynamics  and  thermodynamics  of  the  ocean  are  described 
by  a  set  of  non-linear  basic  governing  equations  containing 
time  and  space  derivatives  requiring  the  definition  of  initial', 
and  boundary  conditions  [Ref. 5].  Usually  we  do  not  have  enough 
information  to  solve  these  conditions  at  smaller  scales,  so  we 
must  choose  a  cut-off  eddy  size  below  which  we  include  only 
the  statistical  effects  of  these  smaller  scales.  The  resulting 
non-linear  equations  will  be  underdetermined  with  no 
analytical  solution.  Therefore,  we  are  limited  to 
approximations  or  numerical  solutions  after  closing  the  system 
with  additional  assumptions  [Ref.2]. 
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1.  Basic  Governing  Equations 

The  governing  equations  are  the  set  of  equations 
describing  the  general  physical  laws  of  fluid  dynamics  and 
their  application  to  the  particular  case  of  ocean  dynamics 
[Ref . 6  and  2 ] : 

a.  Equation  of  State 

This  equation  describes  the  thermodynamic  state  of 
water  relating  the  several  basic  variables  of  density, 
pressure,  temperature  and  salinity.  A  complete  empirical 
relation  in  polynomial  form  is  given,  for  example,  by  Millero 
et.al.  1980  [Ref. 6].  However,  this  equation  is  difficult  to 
handle  and  interpret.  An  approximation  widely  used  is  based  on 
the  fact  that  the  water  is  nearly  incompressible  near  the 
surface  so  that  density  changes  of  a  constant  volume  are  only 
due  to  temperature  or  salinity  changes  (Brian  and  Cox  1972) 
[Ref .5] , 

p  =  p,-y  [l-ax(r-T3)  +px(5-5,)  ] 

where  a  =  2139xlO''’K'^  is  the  thermal  expansion  coefficient  and 
B  =  0.751x10'^  is  the  salinity  expansion  coefficient  and  po  = 
1025.96  kg/m’  is  the  density  correspondent  to  the  temperature 
To  =  15°C  and  salinity  So  =  35  psu  [Ref. 5]. 
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b.  Conservation  of  Mass  (continuity) 


The  conservation  of  mass  for  fluids,  also  called 
the  continuity  equation,  states  that  the  local  rate  of  change 
of  density  is  equal  to  minus  the  mass  divergence,  or: 


dt  dxj 


ipUj)  -0 


If  the  fluctuations  of  density  are  much  smaller  than  the  mean 
density  it  can  be  shown  that  the  incompressibility  assumption 
is  valid  and  this  equation  becomes: 


dp 

dt  dx^ 


=  0 


i.e.  non-divergent  flow  and  constant  density  following  the 
motion . 

c.  Conservation  of  Momentum 

This  balance  describes  Newton's  second  law  of 
motion.  If  we  assume  that  the  vertical  and  horizontal  scales 
are  much  smaller  then  the  earth  radius,  this  conservation  law 
can  be  written  using  the  tangent  plane  approximation  [Ref. 2]: 


dU. 

~dt 


-5,,g-€ 


ep  ^  1 

dx^  p  dx^ 
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where : 


that : 


(5  =  [1  Tn  =  n  €=ri  innq=123  231  312 

-1  m  *  n  1-1  innq=321  213  132 

0  otherwise 

Water  can  be  assumed  to  be  a  newtonian  fluid  so 


"■  . 


dx^  8Xj 


where  /i=l. 4x10-3  kgm'^s"^  and  fit,  are  the  molecular  viscosity  and 
the  bulk  viscosity  respectively  [Ref. 2]. 

However  if  we  assume  that  the  fluid  is  incompressible  and 
defining  kinematic  viscosity  (v)  as  /i/p=l .  4x10*®  m^s'^  we 

obtain : 

1  8x,. 

- _i_  =v - i 

p  dx-  Bx/ 


If  we  assume  that  the  motion  is  shallow,  we  can 
neglect  density  changes  associated  with  the  mass  of  the  fluid, 
though  we  must  retain  their  effect  on  the  fluid  weight,  and  we 
can  further  simplify  the  equations  of  motion  in  the  so  called 
Boussinesq  approximation  [Ref.2]: 


dU, 

it 


dU, 


P 


-c  f  n  - 1 


d^U: 

+  V - i 
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These  equations  represent  the  balance  of  the  inertia  and 
advection  (LHS)  by  gravity,  coriolis  acceleration,  pressure 
gradient  and  viscous  stresses. 

d.  Conservation  of  Heat 

The  conservation  of  heat  describes  the  first  law  of 
thermodynamics.  If  we  neglect  the  change  of  internal  energy 
due  to  viscous  dissipation  and  assume  incompressibility,  it 
can  be  represented  by: 


ar.„  a-'T 

'-fc 


+ 


0 


r 


where  the  left  hand  side  represents  the  storage  and  advection 
of  heat,  Kt=lxl0-7  m^s‘V(pCp)  is  the  coefficient  of  thermal 
diffusivity,  Q  represents  the  rate  of  heat  addition  per  unit 
mass  by  internal  heat  sources  and  Cp=3990  J.Kg’^K'*  is  the 
specific  heat  at  atmospheric  pressure  and  IS^’C  [Ref.  5]. 

The  term  Q  can  be  summarized  by  [Ref. 2, 5  and  6]: 

( 1  -a ) 


where : 


-  short  wave  solar  radiation  (a  =  surface  albedo) 
-  net  heat  gain  by  long  wave  radiation 
Qp  -  net  heat  loss  by  long  wave  radiation 
Qp  -  net  heat  loss/gain  by  conduction 
Q..  -  net  heat  loss/gain  by  evaporation  or  condensation 
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Q^.  -  net  heat  loss/gain  by  advection 
The  terms  Q*  and  will  only  be  significant  at  the  surface 
and  can  be  estimated  by  the  aerodynamical  bulk  formulas,  based 
on  bulk  meteorological  measurements  [Ref.2]. 

The  term  representing  the  heat  associated  with  the  mass 
transport  is  usually  parameterized  by  an  eddy  diffusion 
coefficient  of  the  form: 


0,=K 


>  d-T 

" 


The  terms  associated  with  the  downward  long  wave 
and  short  wave  radiation  are  more  difficult  to  analyze.  We 
must  use  empirical  relations  describing  the  penetration  and 
absorption  of  radiation  with  depth  as  a  function  of 
wavelength . 

Finally  the  heat  loss  by  long  wave  radiation  can  be 
estimated  using  Steffans  law  where  (Qt  «  T*). 
e.  Conservation  of  Salt 

Because  of  the  form  of  the  equation  of  state  used, 
a  salinity  budget  must  be  introduced  of  the  form: 


as  „  ds  d^S 
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stating  that  the  storage  and  advection  of  salt  (LHS)  must  be 
balanced  by  the  molecular  diffusivity  of  salt  and  by  the  net 
salt  gain/loss  by  internal  sources. 

The  term  corresponds  to  the  balance  due  to 
evaporation,  precipitation  and  runoff  effects. 

2 .  Reynolds  Averaging 

The  randomness  found  in  the  smaller  scales  of  motion 
in  a  turbulent  fluid  makes  a  deterministic  approach  to  solving 
these  equations  very  difficult.  If  we  assume  that  the  energy 
at  these  scales  is  the  superposition  of  a  variety  of  small 
size,  band  limited  eddies,  we  can  separate  these  components 
from  the  other  mean  effects  by  averaging  over  a  sufficiently 
large  interval  corresponding  to  a  stationary  low  frequency 
energy  spectral  component.  This  concept  presumes  the  existence 
of  a  spectral  gap  in  the  component  variance  spectra  [Ref.2]. 
However,  the  presence  of  internal  waves,  surface  gravity  waves 
and  many  other  phenomena,  does  not  always  allow  us  to  assume 
the  existence  of  this  gap  between  forcing  and  dissipation 
frequency  ranges.  Nevertheless,  we  can  assume  that  each 
variable  is  composed  by  a  mean  and  a  perturbed  part,  which 
averages  out  to  zero.  However,  the  turbulent  part  of  the 
variance  will  have  non-zero  covariances  with  the  other 
turbulent  components. 
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starting  from  the  previous  basic  governing  eguations 
and  after  separating  the  mean  and  perturbed  components  and 
averaging,  the  resultant  set  of  equations  is: 

1.  equation  of  state 

p  =  p,x[l-ax{T-To)  +px(S-S3)] 


2.  conservation  of  mass  (continuity) 

dV~  _  du'- 
dx.  dx. 


3 .  conservation  of  momentum 


dU,  __  dil. 
dt  ^  dx. 


=  -6^,CT-€,..f  X'.- 


dHJ. 

^ 


a 
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4.  conservation  of  heat 


FL.jtJL.k 

dt  ^  dx,  dx?  F 


(u 


5.  conservation  of  salt 


dS  jj-  -L' 

-r —  C  ,  —z - A  _ - 

dt  '  dx,  '  dx - 
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The  new  terms  introduced  by  this  procedure  are  the 
divergences  of  momentum,  temperature  and  salinity  fluxes,  and 
they  represent  the  interaction  between  the  mean  and  the 
perturbed  flow.  For  these  terms  to  be  significant,  we  must 
have  a  large  Reynolds  number  (  =  advection  terms  / 
molecular  viscosity  terms),  so  that  it  is  usual  to  neglect  the 
viscosity  terms  in  this  set  of  equations.  The  set  of  equations 
represent  an  unresolved  system  of  7  non-linear  equations  with 
22  unknowns  and  are  known  as  the  Reynolds  equations. 

Before  looking  for  any  type  of  solutions  we  first  must 
close  the  system  by  finding  prognostic  equations  for  the  new 
variables  in  the  form  of  fluxes  and  variances,  eventually 
reducing  the  number  of  unknowns.  To  obtain  the  forecast 
equations  for  the  fluxes  and  variances  we  must  first  find 
prognostic  equations  for  the  perturbat j ons .  Following  the 
procedure  described  in  Stull,  1988  [Ref.2],  this  can  be 
achieved  by  subtracting  the  Reynolds  equations  from  the 
governing  equations,  resulting  in  the  system: 


di:  '  :  _  du'  :  ,  U,  ,  du'  ; 
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dr' 

dt  "  dx. 


+  u 


/  dr 

"  dx. 
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dx. 
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In  order  to  obtain  equations  for  the  reynolds  stresses  and  for 
the  scalar  fluxes,  the  following  procedure  was  developed  by 
Stull,  1988  [Ref. 2]: 

1.  Multiply  the  equations  for  a  generic  variable  q' ,  by  u\, 
and  Reynolds'  average. 

2.  Sum  these  resultant  equations,  obtained  by  Reynolds' 
averaging  the  perturbed  momentum  equations  multiplied  by  q,'. 

3.  Use  the  continuity  balance  to  transform  advective  terms 
into  flux  form. 

4.  Neglect  the  coriolis  terms  and  the  pressure,  radiation  and 
molecular  diffusion. 

The  resultant  equations,  are; 
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dY'u  ' ,  _ BT’u' . 


dt 


dx 
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As  we  can  see,  these  equations  do  not  close  the  system  because 
they  introduce  new  variables  in  the  form  of  statistical  third 
moments.  If  we  proceed  looking  for  prognostic  equations  for 
these  third  moments,  they  will  introduce  in  a  similar  way 
fourth  order  type  moments  and  so  forth. 

Therefore  in  order  to  be  able  to  use  these  Reynolds 
equations  additional  assumptions  must  be  used. 

3 .  Turbulent  Closure 

If  there  was  no  turbulence  and  the  statistical  second 
moments  of  the  Reynolds  equations  were  negligible,  they  would 
form  a  closed  system  of  equations  describing  the  mean  motion. 
This  is  typical  of  the  deep  ocean,  where  turbulence  events  are 
episodic  and  of  a  low  mean  level  [Ref. 6]. 
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However,  these  second  moments  are  generally  not 
negligible  on  the  upper  layer  of  the  ocean,  where  atmospheric 
fluxes  can  be  strong.  The  upper  layer,  where  the  Reynolds 
stresses  and  other  turbulent  fluxes  are  important,  is  known  as 
the  turbulent  ocean  mixing  layer.  It  is  not  necessarily  equal 
to  the  well  mixed  upper  layer  of  nearly  constant  density  known 
as  the  remnant  mixed  layer  [Ref. 7].  Therefore,  in  this  study, 
our  closure  problem  is  essentially  a  boundary  layer  problem 
where  mixing  processes  and  air-sea  interactions  take  place. 

In  order  to  close  the  system  of  equations,  we  must 
first  decide  to  what  order  we  want  to  go  to  in  our  prognostic 
equations,  and  at  that  point  parameterize  the  statistical 
moments  based  on  known  variables.  Hence  parameterizations  are 
an  important  tool  that  rely  essentially  on  observations  of 
slow  varying  variables.  However,  a  complete  study  of  these 
parameterizations,  to  see  how  the  observed  variables  relate  to 
the  fast  varying  ones  and  to  validate  the  procedure,  will 
involve  the  simultaneous  observation  at  a  wide  range  of 
scales.  Until  recently  these  type  of  observations  suffered 
from  a  lack  of  instrumentation  capable  of  providing  a  complete 
picture  of  the  ocean  microstructure,  similar  to  those  that 
exist  for  micrometeorology  observations. 

The  turbulent  closure  approximation  may  be  made  to 
order  0,  1,  2,  etc,  depending  on  the  highest  order  moment  in 
the  prognostic  equations.  Some  theories  only  use  a  portion  of 
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the  available  equations  introducing  what  is  known  as  the  half¬ 
order  closures. 

Two  major  approaches  of  turbulent  closures  are  used 
that  are  applicable  at  any  order  [Ref. 2].  In  local  closure,  we 
parameterize  unknown  quantities  at  one  point  by  values  or 
gradients  of  known  quantities  at  the  same  point.  It  assumes 
that  turbulence  has  a  similar  nature  to  molecular  diffusion. 
This  is  like  K-theory,  where  we  assume  the  stresses  are 
directly  proportional  to  the  mean  shear.  In  the  other 
approach,  known  as  nonlocal  closure,  an  unknown  quantity  at  a 
particular  point  is  parameterized  by  known  quantities  at  many 
points  in  space.  It  assumes  turbulence  as  a  superposition  of 
eddies  transporting  fluid  and  its  properties,  like  an 
advective  process.  Examples  of  nonlocal  first  order  closure 
like  the  Transilient  Turbulence  Theory  and  the  Spectral 
Diffusivity  Theory  are  described  in  Stull,  1988  [Ref.2]. 
Usually  higher  order  local  closures  and  nonlocal  closures 
yield  more  accurate  results,  but  they  do  so  by  adding 
computational  expense  and  complexity  [Ref.2]. 

1 .  Local  Zero-Order  Closure 

In  this  approach,  no  prognostic  equations  are 
retained,  not  even  equations  for  mean  quantities.  Mean  flow, 
temperature  and  salinity  are  parameterized  directly  as 
functions  of  space  and  time.  Adimensional  variables  are 
usually  defined  for  certain  classes  of  problems. 
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2.  Local  Half-Order  Closure 


This  closure  scheme  uses  only  a  subset  of  the  first 
moment  equations.  A  variation  of  this  method  yield  the  so 
called  Bulk  Models,  where  temperature  profiles  are  assumed  to 
have  a  certain  shape,  exploring  homogeneity  in  the  mixing 
layer  and  using  vertically  integrated  equations  [Ref.8]. 

3.  Local  First-Order  Closure 

In  this  scheme  we  retain  the  prognostic  equations 
for  the  mean  quantities  and  we  parameterize  the  statistical 
second  moments  (variances  and  cross-variances). 

For  a  generic  variable  ”1"  one  possible 
parameterization  of  its  turbulent  flux  will  be: 


u 


dZ 

dx. 


where  Ki  with  units  m^s'^is  known  as  eddy  viscosity  for  the 
velocity  case  and  as  eddy  diffusivity  coefficient  for  the 
scalar  variables.  This  approximation  is  called  K-Theory  and 
frequently  fails  whenever  large  eddies  are  present  or  when 
there  is  not  a  dominant  scale  [Ref. 2].  A  useful  interpretation 
of  these  approaches  can  be  made  by  looking  at  analogies  to  the 
viscous  stress  expressions  for  Newtonian  fluids.  We  can 
express  the  Reynolds  stresses  in  terms  of  the  shear,  by 
substituting  the  molecular  viscosity  coefficient  with  an  eddy 
viscosity  coefficient,  obtaining: 
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^  reynolds 


=  pK, 


dU 


Since  turbulent  mixing  is  much  more  effective  than 
viscosity,  the  eddy  viscosity  should  be  much  greater  then  the 
molecular  viscosity  [Ref.2].  Also  while  molecular  viscosity  is 
a  property  of  the  fluid,  the  eddy  viscosity  is  a  property  of 
the  flow  so  that  this  coefficient  should  vary  as  turbulence 
level  varies.  Good  parameterizations  should  rely  on  the  Monin- 
Obukhov  length  scale,  on  the  Richardson  number  (Ri)  and  on  the 
stratification,  as  we  can  see  in  the  examples  presented  in 
Stull,  1988  [Ref.2]. 


Another  useful  interpretation  of  this  theory  is  the 
so  called  Mixing-Length  Theory  derived  by  Prandtl  (1925) 
[Ref.2]  for  neutral  conditions.  If  a  turbulent  eddy  moves  a 
water  parcel  upward  by  a  distance  z'  and  if  we  assume  that 
there  is  no  other  mixing  or  changes  on  "1"  within  the  parcel, 
then  the  perturbation  of  ”1"  can  be  written  as: 


J'=J 


SiirioUTidi 


1 


parcel 


similarly 
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In  order  to  have  upward  motion  we  must  have  some  vertical 
velocity  w'.  If  turbulence  is  such  that  it  can  be  assumed  that 


c.  u 


w'--c .  u' 


for 

for 


dz 

ffu 

dz 


<0 

>0 


where  "c"  represents  a  positive  constant.  Then  we  can  define 
mixing  length  "Ml"  as 

Ml  '^-cz'^ 


A  vertical  eddy  flux  of  the  scalar  variable  "1"  will  have  the 
form: 


cz'  dz 


so  that  iC,=AfJ^!4^i 

dz 


Therefore  Ki  should  increase  as  shear  increases  or  turbulence 
increases  and  as  the  ability  of  turbulence  to  cause  mixing 
increases.  Finally  a  remark  should  be  made  about  the 
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assumption  of  linear  gradients  valid  only  over  small 
distances,  noting  that  this  theory  is  a  small-eddy  theory. 

4.  Turbulence  Kinetic  Energy  (TKE) 

TKE  measures  the  intensity  of  turbulence  and  is 
directly  related  to  the  momentum,  heat  and  salt  transport 
through  the  mixing  layer.  Its  balance  is  important  in  the 
physical  interpretation  of  the  fluxes  that  are  present  on  the 
Reynolds  equations,  and  in  describing  the  physical  processes 
that  generate  and  maintain  turbulence.  The  total  kinetic 
energy  per  unit  mass  (KE)  is  given  by: 


using  Reynolds  averaging  we  can  identify  the  following 


components  on  the  mean  KE: 


K=Y'  —1''-'  KE  of  Che  mean  flow 
;  =1 


—  u'l  KE  of  turbulence  {TKE) 
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From  the  Reynolds  equations  we  can  obtain  the  budget  for  the 
KE  of  the  mean  flow  on  the  form  [Ref.2]: 
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Interesting  features  in  this  equation  are  the  fact  that 
Coriolis  has  no  effect  on  the  TKE,  showing  that  it  is  only  a 
redistribution  term.  Another  important  feature  is  the  explicit 
form  of  the  Reynolds  stresses  representing  the  transformation 
of  energy  between  the  large  scales  of  the  mean  flow  and  the 
smaller  scales  of  turbulence. 

From  the  definition  of  potential  energy  (tp)  and  by 
using  the  Reynolds  averaging  procedure  we  can  obtain  the 
following  potential  energy  budget: 


dt 


'q 

P 


We  can  see  from  this  budget  that  the  scalar  turbulent  fluxes 
represent  the  transformation  between  the  mean  potential  energy 
and  turbulence. 

After  identifying  the  role  of  the  Reynolds  stresses 
and  scalar  fluxes  it  is  important  to  see  how  they  balance, 
interact  and  develop  in  time.  This  can  be  accomplished  by 
establishing  a  TKE  budget.  Using  the  prognostic  equations  for 
the  perturbations  and  using  an  analogous  procedure  to  that 
wich  we  used  to  obtain  the  prognostic  equations  for  the 
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Reynolds  stresses  and  scalar  fluxes,  we  can  obtain  prognostic 
equations  for  the  variances  of  the  perturbed  velocity 
components  that  yield  the  following  TKE  balance  [Ref. 2]: 


dt 
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The  LHS  of  the  budget  represents  the  storage  and  advection  of 
TKE  by  the  mean  flow.  This  latter  term  is  often  assumed 
negligible  by  assuming  homogeneous  turbulence  over  the  mixing 
layer . 

The  first  term  on  the  RHS  represents  the  buoyancy 
production/consumption.  Corresponds  to  the  transformation  of 
TKE  into  potential  energy  and  vice  versa.  Usually  it  is 
positive  during  the  day  and  negative  during  the  night  due  to 
diurnal  solar  heating. 

The  second  term  represents  the  mechanical  or  shear 
production  and  is  directly  related  to  the  Reynolds  stresses, 
reflecting  any  energy  sources  due  to  stirring  or  mixing  and 
the  interaction  with  the  kinetic  energy  of  the  mean  flow. 

The  third  term  represents  the  eddy  transport  of  TKE, 
acting  as  a  source  or  loss  whenever  there  is  a  convergence  or 
divergence . 

The  fourth  term  represents  the  pressure  correlation 
term,  describing  how  TKE  is  redistributed  by  pressure 
perturbations.  It  is  this  term  that  is  responsible  for  the 
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interaction  of  turbulence  with  internal  or  surface  gravity 
waves.  The  small  scale  perturbations  can  be  due  either  to 
waves  or  turbulence  and  from  linear  wave  theory  (u'^p') 
represents  energy  flux  or  power;  therefore  we  can  expect 
interaction  between  these  two  phenomena  to  exist  by  sourcing, 
sinking  and  redistributing  TKE. 

The  last  term  represents  the  viscous  dissipation  rate 
or  the  transformation  of  TKE  into  heat. 

B.  MIXING  LAYER  ANALYSIS 

Mixed  layer  depth  can  be  defined  as  the  depth  at  which 
density  increases  by  some  predetermined  amount  above  the 
surface  value,  while  mixing  layer  depth  can  be  defined  as  the 
depth  at  which  the  dissipation  rate  drops  monotonically  below 
some  arbitrary  small  value  [Ref.S].  This  definition  for  mixed 
layer  is  useful  for  observational  purposes,  but  do  not  reflect 
the  depth  range  of  mixing,  i.e.  the  layer  that  is  dynamically 
unstable.  Static  stability  is  a  measure  of  the  capacity  for 
buoyant  convection  that  is  independent  on  the  flow  and  based 
only  on  density  profiles.  Dynamic  stability  measures  the 
tendency  for  mean  flow  to  break  into  smaller  scales,  even 
under  statically  stable  conditions  [Ref. 6].  While  statically 
unstable  flows  become  or  remain  turbulent,  statically  stable 
flows  may  be  laminar  or  turbulent  depending  on  many  factors 
that  can  be  interpreted  based  on  the  TKE  balance. 
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We  can  view  the  dynamic  stability  in  terms  of  a  balance 
between  the  two  different  effects  of  stratification  and  shear 
whenever  the  inertia  is  much  greater  then  the  viscous  effects 
(high  Rg)  and  the  reynolds  stresses  and  scalar  fluxes  are 
significant.  While  stratification  will  tend  to  damp  out  any 
perturbation,  the  shear  of  the  mean  flow  will  tend  to  enhance 
perturbations  by  supplying  the  required  energy. 

A  description  of  a  mechanism  for  the  dynamic  generation  of 
turbulence  is  summarized  by  Stull,  1988  [Ref. 2]  based  on 
laboratory  experiments.  He  describes  the  following  sequence  of 
events : 

1.  shear  across  a  density  interface  (laminar  flow) 

2.  a  critical  value  of  shear  is  reached  and  flow  becomes 
dynamically  unstable,  generating  gentle  waves  at  the  interface 

3.  waves  grow  in  amplitude  and  eventually  break 

4.  between  waves  crests  lighter  fluid  entrains  the  denser 
one  and  patches  of  static  instability  results 

5.  this  static  instability  combines  with  continuous 
dynamical  instability  causing  each  wave  to  become  turbulent 

6.  turbulence  at  the  interface  causes  mixing  of  the  two 
different  fluids  (entrainment) 

7.  mixing  can  reduce  shear  below  critical  value  and 
eliminate  dynamical  instability 

8.  if  there  is  no  continuous  supply  of  energy  to  restore 
the  shear,  turbulence  will  decay  and  flow  will  become  laminar 
again . 
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This  conceptual  model  for  generation  of  turbulence  seems 
consistent,  however  due  to  the  episodic  and  3  dimensional 
nature  of  turbulent  events  in  the  ocean,  observations  of  all 
the  processes  are  very  difficult  to  perform.  Also  the  lack  of 
microstructure  instrumentation  able  to  directly  observe  the 
velocity  field  is  a  major  limitation  on  our  understanding  of 
these  processes. 

Usually  we  rely  on  nondimensional  parameters  to  compare 
these  factors,  like  the  Reynolds  number  and  the  Richardson 
numbers.  The  Reynolds  number  compares  the  inertial  and  the 
viscous  terms  in  the  equations  of  motion,  evaluated  when  the 
Reynolds  stresses  are  more  important  then  viscous  stresses. 
This  happens  only  when  the  non-linear  effects  are  significant. 
For  turbulence  to  exist,  it  is  therefore  necessary  to  have  a 
high  Reynolds  number,  but  this  by  itself  it  is  not  a 
sufficient  condition,  because  it  makes  no  requirements  about 
the  supply  of  energy  necessary  for  turbulence  to  persist. 

The  Richardson  number  compares  the  buoyancy  production  to 
the  mechanical  production  terms  of  the  TKE  balance,  and 
evaluates  when  the  shear  will  be  able  to  overcome 
stratification  to  allow  dynamical  instabilities  to  exist 
[Ref.2].  Usually  we  have  high  Reynolds  numbers  in  the  ocean; 
therefore,  the  Richardson  number  is  the  critical  criteria  to 
determine  if  the  flow  is  turbulent  or  not. 

The  flux  Richardson  number  is  defined  as  [Ref.2]: 
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Rr=- 


_g 

X 


^  ( w^) 


^  dx^ 


For  Rf  less  then  zero  we  have  static  instability  (hence 
turbulent)  flows.  R,  equal  to  zero  corresponds  to  neutral 
flows.  The  critical  conditions  below  which  flows  are  turbulent 
is  Rt  =  1 . 

Using  K-theory,  we  can  parameterize  the  fluxes  based  on 
the  mean  shear  and  obtain  the  gradient  Richardson  number 
[Ref . 2 ] : 


R,= 


p  dz 


The  critical  value  below  which  laminar  flows  become  turbulent 
is  accepted  as  0.25,  while  the  critical  value  for  turbulent 
flows  to  become  laminar  is  1  [Ref.2].  This  hysteresis  is  due 
to  the  fact  that,  for  turbulence  to  start,  we  need  the 
instability  and  a  trigger  mechanism;  while  instabilities  will 
start  at  Rj  near  1,  the  trigger  mechanism  will  only  occur  near 
R,  =  0.25  [Ref .2]  . 
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Because  usually  we  are  limited  to  discrete  observations  at 
a  certain  depth,  the  bulk  Richardson  number  is  defined  as 
[Ref . 2 ] : 

^  p[Alr-Av^] 


This  discretization  will  change  the  critical  values  but  the 
same  critical  value  of  0.25  is  generally  used. 

From  this  discussion  we  can  conclude  that  direct 
microstructure  observations  are  very  important  on  the 
characterization  of  the  mixing  layer.  Mixing  processes  have  a 
special  importance  in  coastal  processes,  controlling  the  rate 
of  cooling  or  warming  of  the  surface  waters  and  the  dispersion 
or  diffusion  of  other  scalar  quantities.  The  role  of  vertical 
turbulent  mixing  on  the  coastal  upwelling  systems  is 
determinant  in  the  characterization  and  understanding  of  these 
type  of  phenomena,  as  we  can  see  in  the  paper  by  Brink  (1980) 
[Ref.9].  Therefore  mixing  layer  observations,  are  need  in 
order  to  study:  mixing  layer  deepening  or  entrainment;  and  air 
sea  interactions,  as  we  can  see  in  the  papers  by  Dewey  et.al. 
and  Mourn  et.al.  [Ref. 7  and  10],  characterizing  coastal  eddies 
and  the  enhancement  of  ocean  fronts  by  turbulent  mixing,  and 
in  the  paper  by  Uwe  et.al.,  analyzing  upwelling  relaxations 
during  the  CODE  experiment  [Ref. 11]. 
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C.  SPECTRAL  DESCRIPTION  OF  TURBULENCE 

The  spectrum,  defined  us  the  Fourier  transform  of  the 
autocorrelation  function,  gives  a  picture  of  how  variance  or 
energy  is  distributed  over  frequency  or  wavelength.  An  eddy  of 
a  certain  wavenumber  "k"  can  be  defined  as  a  linear 
superposition  of  the  disturbances  containing  energy  in  the 
vicinity  of  "k"  [Ref.l].  These  eddies  will  be  exposed  to  the 
strain-rate  field  of  larger  ones,  and  nonlinear  interactions 
and  consequent  energy  transfers  are  expected.  As  a  resul'*',  the 
eddies  will  loose  their  identity  within  one  or  two  periods  or 
wave" jngths . 

The  Fourier  Transform  is  a  decomposition  into  wave 
components  of  different  wavelengths  associated  with  the 
fourier  coefficients.  An  eddy,  however,  has  many  fourier 
coefficients  changing  in  time.  Therefore  more  sophisticated 
transforms  and  a  nonstationary  approach  should  be  used  if  we 
want  to  decompose  the  velocity  field  into  individual  eddies. 

When  analyzing  the  wavenumber  spectra,  i.e.  how  energy  or 
variance  is  distributed  in  scale,  we  can  expect  to  find 
contributions  down  to  the  sub-centimeter  scales  of 
dissipation,  after  which  the  mechanical  energy  is  transformed 
into  heat  by  viscous  processes,  and  the  energy  drops  down  to 
zero  [Ref.l].  If  our  data  is  not  sampled  below  these  scales 
the  spectra  will  be  aliaised.  Therefore  to  understand  both 
momentum  fluxes  and  dissipation  we  must  measure  a  wide  range 
of  scales,  including  sub-centimeter  dissipation  scales. 
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a.  Velocity  Spectrum 

Large  eddies  are  affected  by  the  mean  flow,  such 
that  they  have  a  steady  anisotropy,  whereas  no  permanent 
anisotropy  is  induced  at  an  appropriately  small  scales.  As  the 
eddy  size  becomes  smaller,  we  can  expect  that  the  strain-rate 
field  becomes  isotropic  in  a  mean  sense.  This  concept  is 
called  local  isotropy  and  it  can  not  exist  if  the  Reynolds 
number  is  not  large  enough  (at  least  100  [Ref.l]).  When  local 
isotropy  prevails,  the  time  scales  are  much  shorter  then  those 
of  the  mean  flow,  so  that  small  eddies  respond  quickly  to 
changing  conditions.  This  area  of  the  spectrum  is  called  the 
equilibrium  range  and  includes  all  high  wavenumbers.  In  this 
area  of  the  spectra  the  scales  are  too  small  to  receive  energy 
directly  from  the  mean  flow.  Because  energy  is  finally 
dissipated,  the  amount  of  energy  arriving  to  this  scale  should 
correspond  to  the  total  energy  being  dissipated  by  viscosity. 
The  smallest  scale  that  can  be  reached  on  this  range  is  called 
the  Kolmogorov  microscale  defined  as  [Ref.l]: 


where  v  is  the  kinematic  viscosity  coefficient  and  e  is  the 
dissipation  rate  given  by  [Ref.l]: 


30 


e  =  2.vfk^.S.dk 
0 


and  "S”  is  the  observed  velocity  spectrum. 

For  small  wavenumbers,  viscosity  is  not  important. 
The  important  parameters  to  be  considered  are  those  relating 
transfer  of  energy  from  the  mean  flow  to  turbulence,  and  the 
energy  transfer  to  the  smaller  scales.  The  area  of  the 
spectrum  corresponding  to  the  superposition  of  theL,e  scales 
with  the  Kolmogorov  spectrum  is  often  called  the  Inertial 
Subrange  and  exists  only  for  large  Reynolds  numbers  (bigger 
then  10®).  At  this  range,  no  energy  is  added  or  dissipated,  so 
that  the  energy  flux  across  each  wavenumber  is  constant.  It  is 
characterized  by  a  slope  of  -5/3  on  a  loglog  plot  and  is 
scaled  by  dissipation.  Therefore,  from  observed  spectral 
estimates  in  this  range,  we  can  estimate  the  total 
dissipation . 

Spectral  estimates  will  depend  only  on  the 
wavenumber  and  on  the  total  dissipation.  Using  similarity 
theory  [Ref.2],  we  can  expect  the  spectral  estimate  to  be: 

S(k)  =a.j,.e^^Kk--'^ 

[5]  =Tn^s 
where  [e] 

[ic] 
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where  a„  is  known  as  the  Kolmogorov  constant  with  ranges  from 
1.53  to  1.68  [Ref. 2,  pp.390]. 

For  time  spectra  of  velocity,  two  cases  should  be 
considered.  The  Lagrangian  time  spectra  related  to  the 
temporal  evolution  seen  by  an  observer  moving  with  the 
turbulent  velocity  fluctuations  and  the  Eulerian  spectrum 
related  to  the  evolution  at  a  fixed  point  where  the  mean 
velocity  is  zero.  Their  shape  is  usually  different  due  to 
advection,  but  we  can  still  define  the  inertial  and 
equilibrium  sub-ranges  in  a  similar  way. 
b .  Temperature  Spectrum 

The  temperature  spectrum  is  defined  in  a  similar 
way  to  the  velocity  spectrum,  and  shows  the  same  nonstationary 
character,  with  energy  being  added,  dissipated  and  transferred 
between  different  scales.  If  the  Reynolds  number  is  high 
enough,  there  is  also  an  equilibrium  range,  an  inertial 
subrange  and  an  inertial-convective  subrange  where  the 
temperature  fluctuations  are  simply  convected  and  where  the 
spectral  transfer  of  temperature  variance  is  constant.  In 
certain  cases  when  thermal  diffusivity  becomes  important 
within  the  inertial  sub-range  and  where  viscosity  is  not  yet 
important,  there  is  a  Inertial-Diffusive  subrange;  this  occurs 
in  fluids  with  small  Prandtl  number  (ratio  of  kinematic 
viscosity  to  heat  conductivity).  For  large  Prandtl  number 
fluids,  like  water,  viscous  effects  occur  first  and  then 
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thermal  diffusive  effects;  therefore  we  can  also  define 
Viscous-Convective  and  Viscous-Diffusive  subranges, 
c.  Batchelor  Spectrum 

The  Batchelor  spectrum  corresponds  to  a  theorical 
shape  of  the  wavenumber  spectrum  of  the  temperature  gradient 
fluctuations  for  homogeneous,  locally  isotropic  turbulence  in 
a  homogeneous  fluid  with  a  large  Prandtl  number  [Ref. 12  and 
1].  It  is  defined  by  the  temperature  variance  dissipation 
rate,  molecular  diffusivity  and  by  the  Batchelor  wavenumber 
defined  as: 

_i 


This  wavenumber  as  been  identified  as  the  cutoff  wavenumber 
[Ref . 12  ] . 

The  effect  of  stratification  was  studied  by  Dillon 
and  Caldwell , 1980  [Ref. 12],  who  proposed  that  its  effect 
should  be  more  noticeable  in  the  low  viscous-convective 
wavenumber  ranges  when  turbulence  is  weak  compared  with 
stratification.  One  of  the  most  important  applications  of  this 
theorical  spectrum  is  the  estimation  of  the  kinetic  energy 
dissipation  rates  from  the  spectral  cutoff  wavenumber,  as  we 
can  find  in  the  work  of  Dillon  and  Caldwell , 1980  [Ref. 12]. 
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III.  DOPPLER  SPECTRUM  ESTIMATION 


The  bistatic  velocity  estimation  that  will  be  discussed  in 
the  next  chapter,  requires  a  fast  but  accurate  doppler 
spectral  estimator.  This  estimator  must  rely  on  a  small  number 
of  data  points  and  be  reliable  even  for  weak  signals  submerged 
in  strong  noise.  Several  spectral  estimators  will  be  analyzed 
and  Monte-Carlo  simulations  performed  in  order  to  study  their 
performances . 

Spectral  estimation  corresponds  to  the  estimation  of  a 
function  that  is,  by  definition,  the  Fourier  Transform  of  the 
autocorrelation  function.  The  performance  of  a  spectral 
estimator  should  be  based  on  how  well  we  are  estimating  this 
function  over  a  continuous  range  of  frequencies  [Ref. 13]. 

The  use  of  Spectral  Estimators  to  detect  the  presence  of 
signals  superimposed  in  white  random  noise  demands  a  trade  off 
between  detectability  and  confidence  of  the  detected  values. 
Usually  confidence  is  obtained  by  dividing  the  whole  data  set 
into  smaller  subsets,  estimating  the  spectrum  of  each  subset 
and  then  averaging  the  ensemble  of  estimates.  If  we  consider 
that  each  spectral  estimate  corresponds  to  2  Fourier 
coefficients  and  therefore  has  2  degrees  of  freedom,  the 
ensemble  average  based  on  N  subsets  will  increase  the  number 
of  degrees  of  freedom  to  2N,  improving  the  viability  of  the 
overall  estimate.  However  the  spectrum  based  on  these  subsets 
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will  have  less  resolution  and  the  spectral  estimate  will  be  a 
smoothed  but  more  accurate  version  of  the  unsegmented  one 
[Ref .14] . 

There  are  applications  where  we  need  only  the  behavior  of 
the  spectral  estimator  at  a  specific  frequency  or  in  a 
specific  region,  such  as  the  frequency  estimation  of  a 
sinusoidal  signal  superimposed  in  white  noise.  In  this  case, 
the  smoothed  version  of  the  spectrum  will  be  a  poor  estimator 
of  the  sinusoidal  frequency  as  we  are  seeking  a  parameter 
estimation  (and  not  spectral  estimation),  that  involves 
different  statistics  [Ref. 13]. 

The  frequency  estimator  based  on  a  finite  number  of  data 
points  represents  an  estimate  of  the  true  frequency  that  we 
would  have  if  we  had  an  infinite  number  of  points.  Therefore 
we  need  to  characterize,  in  statistical  terms,  the  behavior  of 
these  estimates.  In  this  work,  the  Bias  and  Variance  of  the 
frequency  estimates  will  be  the  statistics  that  we  will  use  to 
compare  the  performance  of  several  estimators. 

Because  determination  of  these  statistics  for  a  wide 
variety  of  algorithms  is  not  analytically  treatable,  Monte 
Carlo  simulations  are  used  to  obtain  estimates  of  these 
statistics  in  terms  of  the  Signal  to  Noise  Ratio  (SNR)  of  the 
input  signal  and  of  the  record  length  used  in  the  estimator. 
The  results  of  this  simulation  will  be  used  to  determine  the 
appropriate  doppler  frequency  estimator,  for  the  bistatic  data 
stream. 


35 


A.  ESTIMATION  OF  A  SINGLE  COMPONENT  COMPLEX  SIGNAL 
SUPERIMPOSED  ON  GAUSSIAN  WHITE  NOISE 
1 .  Peak  FFT 

This  approach  can  be  interpreted  as  an  estimation  of 
the  parameters  of  a  single  component  sinusoidal  signal  of 
frequency  f ,  superimposed  in  gaussian  white  noise,  of  the  form 
[Ref .15] : 

x(n)  =A  EXP(j2nfn  +  <ff)  *z{n) 

where  z(n)  corresponds  to  complex  gaussian  white  noise  and 
and  "f”  (amplitude  , phase  and  frequency),  assumed 
constant  but  unknown,  are  the  parameters  to  be  estimated. 

The  probability  distribution  function  (PDF)  of  x(n) 
corresponds  to  a  shift  in  the  multivariate  gaussian 
distribution  function  associated  with  the  white  noise,  and  has 
the  form: 


n'^'det 

where  "R„"  is  the  noise  correlation  matrix  and  "s"  the 
sinusoidal  signal.  Since  we  are  assuming  white  noise,  "R„” 
will  be  a  diagonal  matrix.  The  sinusoidal  signal  can  be 
represented  by  the  vector: 

R=A  EXP(j^)  [1  EXP(j2Jlf)  EXPijZUfZ)  .  .  .  EXPijZUfiN-l)]^  = 
EXP[j^)  £ 
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where  e  is  a  vector  containing  the  values  of  the  complex 
exponential  of  frequency  f,  at  the  several  values  of  n. 


To  obtain  a  maximum  likelihood  estimator  (MLE) ,  we 
need  to  find  the  maximimum  of  the  PDF,  which  will  be 
equivalent  to  minimize  the  hermitian  form: 

5=  =  (X-A  EXP(j^)£)^'i2i-A  EXP(j<t>)£) 

First,  if  we  minimize  in  terms  of  the  complex  amplitude 
A^.=A  EXP(j(p),  we  obtain  [Ref.  15]: 

A.  =  ^  x{n)  EXF{-j2Ufr.) 


which  corresponds  exactly  to  the  discrete  fourier  transform  of 
the  data  evaluated  at  f  (assumed  known  in  this  case).  Now  if 
we  substitute  this  expression  into  the  hermitian  form,  we 
obtain : 

N 

to  minimize  this  expression  in  terms  of  the  frequency  we  need 
to  maximize: 
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s-i 


=  —  ly^  x(n)  EXP{-j2'n.fn)  P 


Therefore  the  MLE  of  the  frequency  will  be  the  frequency  where 
the  square  of  the  absolute  value  of  the  discrete  fourier 
transform  (periodogram)  attains  its  maximum. 

In  this  work,  we  used  a  128  point  FFT  algorithm  to 
estimate  the  periodogram  where  the  peak  frequency  will 
correspond  to  the  MLE  estimator. 

This  procedure  will  limit  our  resolution  to  7.8x10'^ 
of  the  sampling  frequency,  so  an  interpolating  procedure 
should  be  used  in  order  to  increase  the  estimator  resolution. 
Nevertheless  this  estimator  suffers  from  the  fact  that,  for  a 
small  number  of  points,  discretization  errors  are  highly 
emphasized  and  its  performance  degradated  if  one  wants  to 
obtain  faster  estimations. 

2.  Pulse  Pair 

Traditionally,  the  Pulse  Paj.r  algorithm  is  used  to 
estimate  the  mean  frequency  of  a  time  series  with  a  single 
frequency  component  signal,  which  corresponds  to  the  frequency 
of  maximum  energy  [Ref. 16].  It  is  based  on  a  very  simple  and 
fast  algorithm  that  has  given  accurate  results  even  for  data 
sets  with  a  small  number  of  points.  However,  white  noise 
contamination  will  induce  a  deviation  from  the  peak  value 
toward  the  zero  frequency  as  the  SNR  decreases. 
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The  algorithm  uses  the  direct  evaluation  of  the  phase 
from  the  autocorrelation  function  estimate.  In  this  project 
only  the  phase  of  the  first  lag  autocorrelation  function  was 
used  and  computed  using  a  classical  approach  based  on 
calculating  covariances  of  two  consecutive  points  averaged 
over  the  sequence  length.  Other  valid  correlation  estimate 
methods  could  be  used  such  as  the  Modified  Covariance  Method 
[Ref. 15],  which  has  the  property  of  giving  estimates  without 
running  out  of  the  data  set.  It  will  be  shown  later,  that  this 
method  is  equivalent  to  an  eigenvalue  decomposition  of  an 
order  2  correlation  matrix  (MUSIC  -order  2  fast  algorithm). 

Higher  lag  autocorrelation  estimates  can  be  used 
either  to  improve  the  estimates  when  the  frequencies  are  much 
smaller  then  the  sampling  frequency,  or  to  obtain  more 
credible  estimates.  However  this  procedure  may  create  problems 
with  the  uniqueness  of  the  estimate,  and  will  be  limited  by  a 
lag,  corresponding  to  the  Nyquist  rate  of  the  highest 
frequency  present  on  the  detected  signal,  after  which  aliaised 
frequency  estimates  occur  [Ref. 16]. 

3 .  Zero  Cross 

This  method  consists  of  the  evaluation  of  a  dominant 
frequency  by  the  number  of  zero  crosses  or  the  number  of  sign 
changes  on  a  zero  mean  time  series.  When  working  with  complex 
signals,  we  must  start  by  defining  criteria  for  a  zero  cross. 
A  point  in  the  real  domain  corresponds  to  a  line  parallel  to 
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the  imaginary  axis  in  the  complex  plane.  Then  using  this 
generalization,  we  can  define  a  zero  cross  (real)  as  the 
change  of  sign  on  the  real  components  of  the  time  series.  On 
other  hand,  a  point  in  the  imaginary  domain  is  a  line  parallel 
to  the  real  axis  in  the  complex  domain;  then  a  zero  cross 
(imaginary)  corresponds  to  a  change  in  sign  in  the  imaginary 
components  of  the  time  series.  Therefore,  we  can  define  a 
complex  zero  cross  as  a  change  in  sign  either  in  the  real  or 
imaginary  part  between  consecutive  points;  therefore  the  zero 
crosses  can  either  be  determined  in  the  real  or  imaginary 
parts.  To  get  a  better  estimate  we  can  also  average  both 
estimates . 

By  itself,  this  simple  method  is  highly  sensitive  to 
high  frequency  noise  contaminations  that  will  increase  the 
number  of  zero  crosses.  This  effect  can  be  partially 
compensated  by  choosing  a  hysteresis  after  each  cross  in  order 
to  minimize  the  high  frequency  contaminations  [Ref. 17].  In 
order  to  apply  this  technique,  we  must  know  the  maximum 
possible  frequency  that  we  can  have  in  the  signal  after  which 
no  frequency  will  be  detectable. 

Another  approach  developed  by  B.Kendem  [Ref. 18] 
suggests  the  incorporation  of  low-pass  filters  and  periodogram 
analysis  to  obtain  more  robust  estimates.  This  procedure  was 
not  developed  in  the  present  analysis,  though  some  interesting 
results  may  be  obtained. 
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4.  AR  Modeling  (Burg's  method) 

AutoRegressive  (AR)  modeling  is  one  general  class  of 
parametric  models  that  we  can  use  to  study  time  series  when 
some  previous  information  exists  [Ref. 15].  It  is  based  on  the 
assumption  that  the  value  of  the  record  at  a  particular  time 
is  determined  by  a  linear  combination  of  P  past  values,  where 
P  corresponds  to  the  model  order,  superimposed  in  white  noise 
(u[n]),  in  the  form: 

P 

x[n\  =-^2  x[n-k]  +u[n] 

ic=l 


The  model  coefficients  are  usually  determined  based  on  the 
autocorrelation  function  of  the  series  using  any  of  the 
estimators  (usually  the  modified  covariance  method  gives  the 
best  results)  by  solving  the  YULE-WALKER  equations. 

Another  method  to  compute  the  model  coefficients  is 
based  on  an  extrapolation  of  the  autocorrelation  function  that 
has  maximum  entropy,  so  that  the  effects  of  truncation 
(spectral  smearing)  are  minimized  [Ref. 15].  The  estimated 
coefficients  are  then  computed  based  on  this  assumption 
directly  using  Lagrange  Multipliers.  This  method  is  known  as 
BURG'S  technique  and  is  the  one  that  was  used  in  the  present 
study  [Ref . 15 ] . 

The  spectral  estimator  will  correspond  to  the  z- 
transform  of  x[n]  evaluated  on  the  unit  circle,  or; 
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where  is  the  sampling  frequency,  H(z)  the  z-transform  of 
the  signal  x[n],  A(f)  the  denominator  polinomial  of  H(z)  after 
the  transformation  z=EXP(-j .2.^.f ) ,  a  the  model  coefficients 
in  vector  form,  and  Pw  the  noise  power. 

The  spectral  estimators  will  have  singular  points  for 
a  discrete  type  spectra  (pure  sin’isoids)  corresponding  to 
zeros  on  the  denominator  of  H(z);  therefore,  in  our  case  of 
only  one  pure  complex  sinusoid  in  white  noise,  we  can  just 
look  for  the  frequency  that  will  be  closer  to  a  singularity. 
This  can  be  accomplished  by  looking  for  the  zeros  of  the 
denominator  of  the  general  spectral  estimation  expression  and 
choosing  the  one  that  is  nearest  to  the  unit  circle.  In  this 
form,  the  frequency  estimation  will  be  totally  determined  by 
the  AR  model  coefficients.  This  is  the  method  we  used  for  the 
AR  dominant  frequency  estimation.  Because  of  its  complexity, 
we  should  expect  this  method  to  be  computationally  intensive, 
and  the  results  will  also  depend  on  the  root  finder  that  we 
use. 

5.  MinimuiD  Variance  Method 

This  method  introduced  by  Capon  [Ref.  13]  is  also 
called  the  Maximum  Likelihood  Method.  It  will  describe  the 
distribution  of  variance  over  frequency,  but  the  area  under 
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this  estimate  will  not  corret  -end  to  the  total  variance 


[Ref. 13].  Therefore  it  is  not  a  true  spectral  estimator,  but 
is  particularly  suitable  for  single  frequency  estimation.  The 
algorithm  computes  the  spectral  estimates  in  accordance  to  the 
output  power  of  a  set  of  narrow  band  finite  impulse  response 
(FIR)  filters  whose  bandwidth  or  resolution  are  dependent  on 
the  frequency  in  an  adaptative  manner  (in  comparison  with  the 
periodogram  that  has  a  fixed  resolution).  The  filter 
coefficients  are  computed  so  that  the  output  variance  is 
minimum  and  its  gain  is  1  (no  distortion).  Ir,  fact  we  do  not 
need  to  compute  all  these  coefficients  explicitly  because  the 
spectral  estimators  can  be  derived  directly  from  the  inverse 
autocorrelation  matrix  in  the  form  [Ref. 13]: 


P 


MT 


_ 1 _ 


where  R  corresponds  to  the  Toeplitz  autocorrelation  matrix, 
that  can  be  estimated  in  several  forms.  In  the  present  work 
the  modified  covariance  method  was  used. 

It  can  be  shown  that  it  gives  spectral  estimates  that, 
for  a  given  order  P,  correspond  to  the  average  of  the  spectral 
estimates  of  the  AR  models  of  order  1  to  P  [Ref.  13],  so  its 
performance  should  result  in  a  smaller  variance.  Its 
resolution  is  between  that  obtained  by  the  classical  spe  tral 
estimators  and  by  an  equivalent  AR  model.  Therefore,  in  order 
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to  obtain  the  single  frequency  estimate,  we  just  need  to 
compute  the  maximum  of  the  above  expression.  However,  due  to 
its  complexity  and  dependence  on  the  resolution  used  in  the 
maximum  value  determination,  we  should  expect  this  algorithm 
to  be  even  more  computational  intensive  than  an  equivalent 
order  AR  method. 

6.  Multiple  Signal  Classification  (MUSIC  Algorithm) 

This  algorithm  is  one  of  several  that  analyze  the 
column  space  of  the  autocorrelation  matrix  into  a  signal 
subspace  and  an  orthogonal  noise  subspace  using  a  Singular 
Value  Decomposition  (SVD)  [Ref. 15].  This  is  accomplished  by 
separating  the  matrix  into  its  eigenvectors  so  that  the 
eigenvectors  with  the  N  largest  eigenvalues,  span  the  signal 
subspace  and  the  others  the  noise  subspace,  where  N  is  the 
number  of  complex  components  of  the  signal. 

This  estimate,  like  the  minimum  variance  method,  is 
not  a  true  spectral  estimator  because  it  does  not  preserve  the 
total  variance,  however  it  will  work  as  a  discrete  frequency 
parameter  estimator.  Its  general  expression  is  [Ref. 15]: 

^MUSIC^  p 
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where  M  is  the  number  of  complex  sinusoids  in  the  signal  and 
v'  are  the  eigenvectors  that  span  the  noise  subspace 
(associated  with  the  smaller  eigenvalues). 

The  discrete  signal  components  will  be  orthogonal  to 
any  linear  combination  of  the  eigenvectors  that  form  the  noise 
base,  and  so  we  could  expect  singularities  at  the  signal 
frequencies.  In  reality,  we  will  not  have  pure  singularities 
but  the  spectral  estimates  will  peak.  For  our  purpose  of 
estimating  a  single  frequency  component,  the  procedure  used 
was  based  on  finding  the  roots  of  the  averaged  noise  space 
eigenvectors  (taken  as  polynomials  of  the  complex 
exponential),  whose  phases  correspond  to  the  peak  frequencies. 
The  eigenvector  associated  with  the  root  closest  to  the  unit 
circle  is  the  single  frequency  estimator. 

It  can  be  proved  that  when  we  use  an  order  2  matrix  in 
a  single  complex  frequency  component  signal,  the  frequency 
estimator  is  equivalent  to  the  Pulse  Pair  algorithm,  and  so 
the  frequency  estimate  can  be  obtained  directly  from  the 
autocorrelation  function. 

The  autocorrelation  function  for  a  single  complex 
component  signal  superimposed  in  white  noise  will  have  the 
form: 

r^ik)  =P.  exp  ( j2nf-^  )  +p„6  (k) 
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and  the  corresponding  2"'^  order  autocorrelation  matrix  will  be: 


The  eigenvector  associated  with  the  smaller  eigenvalue  (-p„) 
will  be: 


i?(o)  y2l^ 


\RiO) 
(1) 


v/2 


exp{-j2n-^) 

^  s 


Therefore,  the  determination  of  the  frequency  estimate  in  this 
case  will  correspond  to  the  determination  of  the  phase 
associated  with  the  first  lag  autocorrelation  function 
estimate,  as  in  the  Pulse  Pair  estimator. 

This  analytical  solution  of  the  2"''  order  eigenvalue 
problem  was  implemented  for  performance  comparisons,  and  a 
singular  value  decomposition  was  used  for  higher  orders.  The 
Modified  Covariance  Method  was  used  to  estimate  the 
autocorrelation  matrix  at  any  order. 

B.  MONTE-CARLO  SIMULATION 

To  evaluate  the  performance  of  these  spectral  estimators 
for  the  estimation  of  the  frequency  of  a  complex  sinusoid 
superimposed  in  complex  independent  gaussian  white  noise,  a 
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MONTE  CARLO  simulation  was  built  within  MATLAB  environment. 
This  simulation  is  based  on  50  independent  trials  where  the 
BIAS,  VARIANCE  and  a  TIME  INDEX  are  evaluated  at  each  step. 
The  model  has  the  input  variables: 

1.  Signal  to  Noise  Ratio  (scalar  or  vector) 

2.  Number  of  Points  (scalar  or  vector) 

3.  Signal  Frequency  (scalar  or  vector) 

4 .  Type  of  Spectral  Estimator 

The  model  output  variables  are: 

1.  Bias 

2.  Variance 

3.  Time  Index  -  proportional  to  the  computation  time 
and  are  stored  at  each  step. 

The  algorithm  of  the  simulation  program  is  as  follows: 

1. SELECT  TECHNIQUE 

1. Peak  FFT 

2.  Pulse  Pair 

3.  Zero  Cross 

4.  Modified  Prony's  Technique 

5. AR  Model  (burg  algorithm) 

6.  MV  (capon's  method) 

7.  Music  (eigenvalue  decomposition! 

2. SELECT  PARAMETERS  f increments . initial  and  final  values) 

1.  Number  points  <->  SNR 

2.  Frequency  <->  SNR 
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3. Frequency  <->  Number  Points 


3. RUN  MODEL 

1.  Initialize  variables 

trials=50 

test  frequency=20  Hz 
sampling  f requency=100  Hz 

2.  Main  loop 

generate  noise 

generate  new  frequency 

generate  signal 

estimate  frequency 

compute  estimate  bias  and  variance 

compute  processing  time 

The  complex  independent  gaussian  white  noise  is  generated 
by  a  function  that  uses  different  seeds  for  the  real  and 
imaginary  parts  at  each  step,  in  order  to  generate  independent 
complex  gaussian  noise.  Each  noise  component  is  then  averaged 
from  5  different  trials  to  obtain  a  pseudo  gaussian  noise. 
This  number  of  5  trials  represents  a  compromise  in  the 
computation  speed.  The  generated  noise  is  normalized  to 
guarantee  a  variance  of  1. 

Model  runs  using  a  20  Hz  signal,  a  sampling  frequency  of 
100  Hz,  over  a  -20  to  10  SNR  and  16  to  128  points  time  series, 
were  made  for  the  following  spectral  estimators: 

1 .  Peak  FFT 
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2.  Pulse  Pair 


3.  Zero  Cross 

4 .  AR  Model  order  4 

5.  MUSIC  order  4 

6.  MUSIC  order  3 

7.  MUSIC  order  2 

8.  MUSIC  order  2/Pulse  Pair  (fast  algorithm) 

The  selected  plots  of  the  bias  and  the  logarithms  of  standard 
deviation  are  shown  on  the  Figures  3.B.1-6. 

1 .  Discussion  of  Some  Results 

The  performance  of  an  estimate  can  be  assumed  good 
when  its  bias  is  less  then  1  Hz  and  the  logic  of  the  standard 
deviation  is  less  then  0.  Using  this  criteria  the  following 
tables  show  the  SNR  at  which  the  bias  of  the  several 
estimators  is  equal  to  1  Hz,  and  the  log  standard  deviation 
equal  to  0: 

TABLE  1  -  SNR  at  BIAS  =  1  Hz 


Estimator 

20  Points 

40  Points 

60  Points 

80  Points 

Peak 

-4 

-5 

-7 

-8 

FFT 

Pulse 

-3 

-8 

-8 

-8 

Pair 

9 


Estimator  I  20  Points  40  Points  60  Points  80  Points 


Zero 

Cross 


AR  Model 
(order  4) 


MUSIC 
(order  4) 


20  Points 

40  Points 

60  Points 

80  Points 

Estimator 


Peak 

FFT 


Pulse 

Pair 


Zero 

Cross 


AR  Model 
(order  4) 


Estimator 

20  Points 

40  Points 

60  Points 

80  Points 

MUSIC 

(order  4) 

10 

5 

2 

2 

MUSIC 

(order  2) 

5 

3 

2 

0 

2.  Conclusions 

For  all  the  estimators,  after  a  sequence  length  (NP) 
between  20  and  64,  the  estimator  variance  becomes  dependent 
only  on  the  SNR.  As  a  general  rule  except,  for  AR  modeling, 
the  variance  is  more  restrictive  then  the  bias.  This  means 
that  we  can  get  some  improvement  in  the  performance  of  this 
estimators  if  we  take  averaged  measurements  of  the  estimated 
frequency.  The  estimator  that  seems  to  work  better  for  low  SNR 
is  the  PEAK  FFT  followed  by  AR  modeling  (order  4),  but  AR 
modeling  has  poor  processing  time  indexes  (too  slow).  The 
MUSIC  estimator  is  very  sensitive  to  the  noise  contamination 
(SNR).  The  equivalent  PULSE  PAIR  method  and  the  MUSIC  (order 
2  fast  algorithm)  do  not  show  any  significant  improvement  by 
changing  the  AUTOCORRELATION  estimators,  at  least  between  the 
two  that  were  used  in  this  program  (biased  and  the  modified 
covariance  estimators  [Ref.l5]). 

The  fastest  algorithm  is  the  MUSIC  (order  2  fast 
algorithm)  at  all  sequences  lengths  followed  by  the  PULSE  PAIR 
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Figurr  l.B.l  -  Montp-CarJo  Simulation:  Peak  FFT  algorithm 
Ptap  and  Variance  plots. 
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Fiqiirc  3.B.5  -  Monte-Carlo  Simulation: 

algorithm,  Bias  and  Variance  plots 
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and  ZERO  CROSS  for  small  sequences  and  the  PEAK  EFT  for  large 
sequences . 

The  lowest  SNR  for  which  we  can  get  unbiased  estimates 
with  a  good  variance  is  -10  when  using  a  PEAK  FFT  with  128 
point  sequences.  Unfortunately  this  performance  indicator  does 
not  show  the  discretezation  error,  which  becomes 
overwhelmingly  significant  for  small  FFT  lengths.  We  can  try 
to  reduce  this  error  by  using  interpolation  schemes,  however 
this  procedure  will  make  this  algorithm  slower  for  a  small 
number  of  points  and  is  very  susceptible  to  low  SNR.  For  SNR 
between  -5  and  0  we  should  use  sequences  of  at  least  32  to  64 
data  points  (PEAK  FFT  also),  but  these  estimates  will  have 
severe  discretization  errors.  For  SNR  between  0  and  5  we  can 
use  either  the  PEAK  FFT  or  the  PULSE  PAIR/MUSIC  2  with  128  to 
32  sequence  length.  For  SNR  greater  then  5  the  best 
performance  is  given  by  the  MUSIC  (order  2  fast  algorithm)  or 
PULSE  PAIR.  If  we  want  to  emphasize  the  speed  of  the 
estimators,  the  best  results  are  the  ones  given  by  the  PULSE 
PAIR/MUSIC  2,  requiring  high  SNR. 

Because  the  CDV  package  bistatic  doppler  frequency 
estimation  is  based  on  32  points  records,  and  the  algorithms 
run  in  real  time  using  C  compiled  code  in  a  DSP,  the  most 
suitable  algorithms  selected  from  these  simulations  are  the 
PEAK  FFT  for  poor  SNR,  and  the  PULSE  PAIR/MUSIC  2  algorithms 
for  higher  SNR.  When  considering  discretization  noise,  the 
first  choice  for  a  spectral  estimator  to  be  used  in  the  CDV 
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real  time  processing  will  be  the  PULSE  PAIR  estimator,  because 
of  its  performance,  simplicity  and  velocity. 
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IV.  THE  CDV  MICROSTRUCTURE  PACKAGE 


A.  RECENT  METHODS  FOR  TURBULENCE  OBSERVATIONS 
1.  Rapid  Sampling  Vertical  Profiler  (RSVP) 

The  RSVP  was  developed  by  Caldwell,  Dillon  and  Mourn  in 
1987  [Ref .19]  for  sampling  near  microstructure  parameters  from 
a  moving  vessel  permitting  the  study  of  the  horizontal 
variability  of  upper  layer  profiles.  The  instrument  makes 
near-microscale  measurements  (3  cm  resolution)  of  temperature, 
salinity  and  turbulence  dissipation  down  to  00  m  when 
deployed  from  a  ship  moving  up  to  6  Knots. 

The  profiler  consists  of  a  small  retrievable  probe  at 
the  end  of  a  cable  that  is  also  a  data  link.  Damping  flaps  are 
used  to  ensure  uniform  vertical  velocity  and  orientation.  The 
airfoil  shear  probe,  used  to  estimate  the  TKE  dissipation 
rate,  is  very  sensitive  to  vibrational  noise  of  the  profiler, 
which  may  dominate  the  signal  in  many  situations.  The  process 
by  which  it  senses  shear  is  based  on  the  deflection  of  a  piezo 
ceramic  beam  in  a  rubber  jacket  by  cross-stream  velocity.  By 
careful  instrument  design,  virtually  all  current  measurements 
of  ocean  turbulent  velocity  have  been  made  with  this  type  of 
sensor . 

Mourn  et.al.  [Ref. 10]  performed  upper  ocean 
observations  with  this  instrument  where  they  were  able  to 
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characterize  the  fine  structure  of  transition  zone  coastal 
eddys,  and  Dewey  et.al.  1990  [Ref.  7]  used  this  type  of 
measurements  to  study  the  enhancement  of  ocean  fronts  by 
turbulent  mixing.  The  interesting  results  obtained  in  both 
experiments  show  the  importance  and  applicability  of  the  RSVP 
in  upper  ocean  observations - 

2.  Pulse  to  Pulse  Coherent  Sonar 

A  pulse  to  pulse  coherent  sonar  was  developed  by 
Lohrmann  et.al.  1990  [Ref. 20],  called  the  High  Resolution 
Current  Profiler.  This  instrument  uses  four  307  KHz 
transceivers,  mounted  30“  off  the  vertical  in  90°  azimuthal 
intervals.  Fifty  range  cells,  covering  a  profiling  distance  of 
10  meters,  measure  radial  current  components  along  the 
transceiver  beams.  These  measurements  do  not  give  the  true 
velocity  vector,  but  assuming  a  stationary  and  homogeneous 
turbulent  field,  allows  an  estimation  of  the  Reynolds  stresses 
and  turbulent  kinetic  energy  by  separating  the  mean  and 
fluctuating  parts  of  the  radial  speeds.  In  this  way,  this 
instrument  remotely  senses  the  turbulent  stresses  without 
inducing  any  flow  distortions  and  provides  moderate  vertical 
resolution . 

3.  Bistatic  Acoustic  Doppler  Profiler 

This  instrument,  developed  by  T.P.  Stanton  (1990), 
measures  the  explicit  velocity  vectors  at  a  very  high  temporal 
and  spatial  sampling  rate  along  one  direction.  The  system 
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consists  of  four,  300  KHz  fan-beam  receiver  transducers, 
radially  separated  from  a  central  narrow  beam-width 
transmitter  transducer,  which  rapidly  emit  short  duration 
acoustic  pulses.  Range  gated  doppler  shifts  are  evaluated  at 
each  receiver  using  dual  interpulse  period,  coherent  doppler 
estimation  methods. 

Simultaneous  receptions  at  the  four  transducers 
provide  an  overdetermined  estimate  of  the  three  velocity 
components  at  range  bins  along  the  ensonified  volume, 
typically  every  10  cm  for  10  m  ranges.  This  system  allows 
concurrent  measurement  of  reynolds  stress  components  and  fine- 
scale  shear  out  to  ranges  of  30  m. 

Velocity  ambiguities  from  the  coherent  velocity 
estimation  limit  the  unambiguous  velocity  range  and  profiling 
range.  Unfortunately,  the  large  frame  for  mounting  the 
acoustic  components  restricts  applications  of  the  system. 

B.  SPECIFICATIONS  OF  THE  CDV  SYSTEM 

The  CDV  Microstructure  Package  is  a  system  developed  by 
T.P.  Stanton  (1991)  to  measure  ocean  turbulence  variables  at 
a  sub-centimeter  resolution.  It  has  3  convergent  acoustic 
transceiver  beams  operating  at  5.296  MHz  in  a  backscatter 
coherent  mode,  using  both  bistatic  (0.6  cm  focal  cube)  and 
monostatic  operating  modes.  In  the  bistatic  mode,  one 
transducer  is  transmitting,  while  the  other  two  are  receiving 
the  scattered  energy  from  a  focal  volume.  In  the  monostatic 
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mode  the  same  transducer  that  transmits  the  signal  receives 
the  backscattered  energy  at  range  bins  along  the  acoustic 
beam . 

The  acoustic  measurements  provide  the  velocity  vector, 
shear  and  kinetic  energy  estimation  at  a  44  Hz  rate,  while 
simultaneous  microstructure  measurements  of  temperature, 
conductivity,  pressure,  shear  and  the  3  components  of  the 
instrument  acceleration  are  made. 

This  package  operates  in  a  profiling  mode  at  approximately 
25  cm/s  descending  or  ascending  rate,  or  fixed  to  a  tripod 
facing  down  for  bottom  layer  observations,  or  facing  up  for 
shallow  water,  near  surface,  observations.  The  package  has 
also  been  attached  to  an  ROV  for  semi-lagrangian  upper  ocean 
turbulence  measurements. 

The  digital  data  is  transmitted  to  the  surface  by  a  1.4 
Mbaud  fiber  optic  data  stream  that  is  stored  on  a  PCM/VCR 
system.  The  data  stream  is  also  passed  to  an  AT&T  digital 
signal  processor  (DSP)  in  a  386  based  HP  VECTRA  for  real  time 
processing  and  graphical  monitoring  of  selected  channels.  This 
pre-processed  data  is  then  stored  on  optical  disks  for  further 
processing.  A  summary  of  the  profiler  package  technical 
specifications  follows: 
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LENGTH 

DIAMETER 

WEIGHT 

1.2  ra 

12  cm 

7  Kg  (in  air) 

TABLE  4  -  ACOUSTIC  CHANNELS 


FREQUENCY 

POWER 

SAMPLING 

OUTPUT 

FREQUENCY 

FREQUENCY 

5.296  MHz 

5  Watt  peak,  in 

11.1  KHz  (after 

44  Hz 

a  2°  beamwidth 

demodulation) 

TABLE  5  -  ANALOG  CHANNELS 


SAMPLING 

TEMPERATURE 

CONDUCTIVITY 

OUTPUT 

FREQUENCY 

PRE-  EMPHASIS 

PRE-EMPHASIS 

FREQUENCY 

1.4  KHZ 

0.339  Hz 

0.100  Hz 

44  Hz 

5,6  KHz  (shear) 

C.  SIGNAL  PRCXTESSING  OVERVIEW 

The  acoustic  inodes  are  controlled  by  a  programmable 
sequencer  that  rotates  the  transmitting  and  receiving  modes  of 
the  three  transceivers,  such  that  3  monostatic  and  3  bistatic 
modes  have  been  selected  every  20  ms.  This  acoustic  data  is 
received,  complex  demodulated  to  baseband  and  the  complex  (in 
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phase  and  quadrature)  channels  sampled  at  11.1  KHz,  for  the 
bistatic  data  and  88.2  KHz  for  each  monostatic  channel. 

The  temperature  and  conductivity  variables  are  pre¬ 
emphasized,  in  order  to  obtain  an  improved  SNR  at  small  scale 
and  directly  estimate  gradient  quantities,  and  sampled  at  1 . 4 
KHz.  The  3  acceleration  channels  and  pressure  are  sampled  at 
the  same  rate,  while  the  shear  pre-emphasized  data  is  sampled 
at  5.6  KHz. 

The  data  stream  comprises  a  set  of  8  word  records,  such 
that  the  first  two  are  used  by  the  non-acoustic  channels  and 
the  other  3  pairs  correspond  to  the  complex  demodulated 
acoustic  data.  These  8  word  records  are  contained  in  a 
sequence  frame  of  2051  words  (corresponding  to  22.7  ms  of 
data)  which  are  passed  by  DMA  to  the  DSP  for  processing  and  to 
a  PCM/VCR  system  for  stor<ige. 

Once  the  DSP  acquires  the  data,  the  data  stream  is 
interpreted  and  the  real  time  processing  is  performed. 
Further  processing,  such  as: 

1  -  Water  salinity  computation 

2  -  Water  density  computation 

3  -  Velocity  vector  (u,  v,  w  components)  estimation 

4  -  Spectral  analysis  of  several  variables 

5  -  Reynolds  stresses  estimation 

6  -  Scalar  fluxes  estimation 

7  -  Dissipation  estimation 

8  -  T.K.E.  estimation 
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is  performed  by  external  programs  accessing  this  preprocessed 
data  on  the  optical  disks.  The  overall  flow  graph  is  shown  in 
the  figure  4.C.I. 

D.  ACOUSTIC  CHANNELS  -  VELOCITY  ESTIMATIONS 

The  velocity  component  estimation  in  the  acoustic  bistatic 
mode  is  based  on  the  mean  doppler  shift  of  the  backscattered 
acoustic  energy  resulting  from  the  relative  velocity  of 
scatterers  present  in  the  focal  volume.  The  real-time  pre- 
processed  acoustic  data  stored  on  the  optical  drive  comprise 
the  mean  doppler  shifts  for  the  6  different  acoustic  paths 
between  pairs  of  transmitters  and  receivers.  In  the  bistatic 
mode,  the  single  transmitter/two  receiver  configuration  is 
sequenced  around  the  three  transducers.  In  order  to  obtain  the 
velocity  components  from  this  spooled  data,  the  following 
steps  are  performed: 

1  -  Clear  data  of  large  spikes  using  a  reciprocal 
channel  comparison  and  a  4’'^  order  linear  predictor,  based  on 
the  previous  20  point  autocorrelation  function  estimate. 

2  -  Compute  velocities  from  the  doppler  frequencies. 

3  -  Transform  estimates  to  an  orthogonal  frame  of  reference 
fixed  to  the  instrument. 

4  -  Rotate  estimate  to  an  absolute  vertical  reference  frame. 

1 .  Doppler  Velocities  Estimation 

For  two  generic  transducers,  where  transducer  1  is 
transmitting  and  transducer  2  is  receiving,  see  figure 
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Fiqurp 


CDV  Signal  Processing  overview. 


4.D.1.1,  and  assuming  that  the  equipment  is  fixed  to  the 
bottom,  we  can  derive  the  following  doppler  relationships 
[Ref. 21].  If  a  pulse  is  transmitted  at  time  t=0,  it  will 
arrive  at  BS  at  time  t=ti,  that  for  a  sound  velocity  ”c"  will 
be  given  by: 

c-Vbs 


Now,  if  a  second  pulse  is  transmitted  at  time  tj,  it  will 
arrive  at  BS  at  time  t,  and 


t3  =  - 


J  +  C.  t. 


c-v. 


T1 


-•tj-  ■ 


BS 


C-V, 


n  •  ^2 


BS 
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Assuming  that  during  the  interval  (Ojtj)  we  are  transmitting 
at  a  frequency  there  will  be  wave  fronts  arriving  at 
BS,  during  the  interval  tj-t^.  Therefore,  the  frequency 
detected  at  BS  will  be: 


^ BS  '  ^ 


=  f  ^  =»f  = 

tx  •  ''2  ^BS 


c-v, 


T1 

BS 


.  f 


Cx 


This  frequency  corresponds  to  the  frequency  of  the  scattered 
signal  at  BS  that  will  be  send  towards  T2 .  Using  the  same 
procedure  we  arrive  at  the  following  expression  for  the 
frequency  at  T2 : 


— £_  .  f  =£_^  .  f. 
C+VfV  C-*-Vps 


If  we  suppose  that  the  velocity  at  BS  has  one  component  V„, 
see  figure  4.D.1.2,  then  =  -Vbs^'^  and  fij  =  ftx*  Therefore 
the  only  velocity  component  producing  a  doppler  shift  will  be 
along  the  direction  defined  by  the  axis  bisection  of  the  two 
beams  i.e.  V^.  This  component  will  have  equal  projections 
along  the  two  beams,  so  that: 


Therefore 
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Fiqure  4.D.1.2  -  Velocity  projections  scheme  for  a  2 

transducers  bistatic  system. 


I’'' =r. 


Af''- _ 


1  r 

2  •  • 


and 


V 


r . 


AX’- 
'  t,: 


Because  of  the  symmetric  geometry  of  this  system,  this 
expression  will  also  be  true  for  the  case  when  T2  is 
transmitting  and  T1  receiving.  This  fact  creates  redundancy  in 
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Figure  4.D.2.1  -  CDV  Package,  Orthogonal  Reference  Frame 

fixed  to  the  instrument. 


the  observed  dopplcL  frequency,  which  provides  increased 
confidence  of  the  estimates  by  averaging,  and  allows  the 
detection  of  single  spurious  wild  points.  Therefore  the  six 
acoustic  transceiver  modes,  will  measure  three  slanted 
velocity  projections  along  the  axes  bisecting  the  acoustic 
beams . 

If  the  instrument  was  in  motion  we  would  have  to  add 
the  instrument  velocity  components  along  these  axis,  in  order 
to  get  absolute  values.  These  operations  are  implemented  by 
programs  external  to  the  real  time  processing. 


2.  Velocities  on  an  Orthogonal  Reference  Frame 

The  Doppler  velocity  components  must  be  transformed 
into  an  orthogonal  frame  before  the  Reynolds  stresses  or  the 
scalar  fluxes  can  be  estimated.  We  can  define  the  instrumental 
orthogonal  reference  frame  as  in  figure  4.D.2.I. 

The  scatterers  velocity  will  be  defined  as: 

y  .=Uji+Vjj*Wjk 


and  the  observed  velocities  will  be: 


cosy 

0 

siny 

cosycosa 

cosycos  (a--^ ) 

siny 

X 

cosycosa 

-cosycos  ) 

siny 

i .  e. 


Therefore,  this  system  may  be  solved  in  terms  of  the  unknowns 
u,  v  and  w  by  inverting  the  non-singular  matrix  "T". 

The  next  step  is  to  rotate  this  velocity  in  the  x-axis 
and  y-axis  in  order  to  obtain  a  vertical  z-axis.  This  can  be 
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accomplished  in  several  ways  depending  on  each  particular 
deployment.  For  the  case  when  we  have  the  system  fixed  to  a 
tripod,  we  can  look  for  the  rotations  that  minimize  the  mean 
vertical  velocity,  or  just  measure  the  platform  orientation 
using  the  3  axis  accelerometer  as  tiltmeters. 

E .  NON- ACOUSTIC  DATA 

The  non-acoustic  channels  on  the  CDV  package  are:  3 
component  acceleration,  pressure,  temperature, conductivity  and 
shear.  The  real  time  processing  for  these  channels  was  coded 
in  C  language  and  compiled  on  the  DSP  together  with  the 
Assembler  routines  that  access  the  CDV  data  stream  in  memory. 
This  DSP  program  is  controlled  by  another  C  routine  running  on 
the  Host  computer,  that  interacts  with  the  user,  reads  the 
proper  setup  files,  receives  the  data  from  the  DSP,  plots 
selected  variables  of  the  pre-processed  data  and  sends  the 
output  to  an  optical  mass  storage  unit.  The  details  on  the 
processing  of  each  particular  channel  are  described  next. 

TEMPERATURE  AND  CONDUCTIVITY 

The  temperature  variable  (T)  is  measured  using  a  high 
frequency  response  FP  07  thermistor,  and  conductivity  (C) 
using  a  2  electrode  AC  conductivity  probe.  The  analog 
temperature  and  conductivity  are  pre-emphasized  at  20  dB  per 
decade  for  frequencies  greater  than  0.339  Hz  and  0.100  Hz 
respectively,  for  SNR  optimization.  They  are  sampled  at  1.4 
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KHz  and  the  multiplexed  data  stream  is  sent  to  the  DSP  by  DMA 
assembly  input  routine. 

A  DSP  function  demultiplexes  variables  from  the  data 
frame  to  separate  selected  variables  from  the  stream.  A  first 
difference  criteria  with  a  threshold  s=t  by  a  setup  file,  is 
used  to  detect  spikes  on  the  data.  Single  bad  points  are 
substiif.uted  by  the  previous  valid  point.  Deterministic 
coherent  noise  over  each  frame  period  (32  words)  is  subtracted 
from  each  data  cycle.  This  coherent  noise  is  computed  from  the 
mean  of  approximately  5000  frames  and  is  entered  from  an  input 
file.  The  raw  T  and  C  gradient  data  are  then  squared  and 
accumulated  on  1  second  rotary  buffers  for  estimation  of 
gradient  variances.  Both  the  T  and  C  data  are  also  de- 
emphasized  by  single  pole  digital  Butterworth  filters  and  low 
pass  (LP)  filtered  by  a  second  order  Butterworth  filter  with 
a  cutoff  frequency  at  20  Hz.  These  operations  are  done  in 
sequence  by  different  functions  to  avoid  round  off  errors. 
Finally  the  de-emphasized  data  is  decimated  to  44  Hz  and  the 
variances  are  computed  based  on  the  mean  of  the  1  second 
rotary  buffers.  These  44  Hz  data  are  sent  every  second  to  a 
HOST  program  running  on  the  HP  Vectra. 


Acceleration  is  observed  using  3  orthogonal  solid 
state  DC  response  accelerometers  in  the  head  of  the  CDV. 
Pressure  is  sensed  using  a  strain  gauge  pressure  sensor.  Both 
Acceleration  and  Pressure  are  processed  in  a  similar  manner. 


The  DSP  function  that  demultiplexes  the  input  buffer,  also 
separates  these  4  channels  from  the  data  stream.  The  same 
deglitching  procedure  is  used  concurrently  with  a  threshold  of 
200  raw  counts  for  the  four  channels.  There  is  negligible 
observed  coherent  noise  on  these  channels.  The  data  is  then  LP 
filtered  by  a  third  order  Butterworth  with  a  cutoff  frequency 
at  20  Hz,  decimated  to  44  Hz  and  sent  to  the  HOST  program 
every  1  second  for  monitoring  and  storage. 

SHEAR 

To  measure  the  velocity  shear,  cross-velocity 
variations  are  measured  with  a  high  sensitivity  airfoil  probe. 
Because  of  its  nature,  this  sensor  is  easily  noise 
contaminated  by  variations  in  the  instrument  velocity  or 
spurious  vibrations.  Also,  it  only  works  when  the  instrument 
is  in  motion  and  senses  only  velocity  changes  perpendicular  to 
the  mean  sensor  velocity. 

The  shear  probe  is  sampled  at  5.6  KHz  and  the  data 
send  to  DSP  by  DMA  assembly  input  routine.  The  same  DSP 
function  separates  the  shear  from  the  data  stream  and  applies 
a  200  count,  first  difference  deglitching  criteria.  Any 
deterministic  coherent,  noise  is  subtracted  from  each  data 
cycle.  Next,  the  pre-emphasized  (gradient)  data  is  squared  and 
accumulated  in  1  second  rotary  buffers,  and  the  variance  is 
computed  at  a  44  Hz  rate  and  sent  to  the  HOST  program. 
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F.  CALIBRATION  AND  ERROR  ANALYSIS 


The  temperature  data  was  previously  calibrated  using  a 
third  order  polynomial  fit  of  the  form: 


TEMP= A.  COUNTS^  *B.  COUNTS'^  *C.  COUNTS 
where 
/9=10 . 662337 
fi=-3. 957935, 10'^ 

C=l. 216833. 10'^ 
£>=-1.345511. 


where  TEMP  represents  temperature  values  in  Degrees  Celsius 
and  COUNTS  the  CDV  outputed  values. 

The  pressure  sensor  was  calibrated  by  recording  the 
pressure  values  obtained  from  water  columns  from  0.10  meter  to 
2.00  meter,  0.10  meter  apart.  The  calibration  transformation 
was  then  computed  using  a  least  squares  fit  to  a  straight  line 
by  solving  a  system  of  20  equations  on  2  unknowns.  The 
resulting  calibration  is: 

DEPTH=0 . 000191964 . COUNTS-0 . 0218414 

where  DEPTH  is  the  water  column  height  in  meters  and  COUNTS 
the  pressure  channel  output  values.  The  mean  square  error  is 
0.0012  m.  The  calibration  was  performed  with  pure  water,  so 
that  in  order  to  compute  the  depth  in  salty  water,  this  values 
must  be  multiplied  by  a  factor  related  to  the  density  anomaly. 
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that  for  practical  use  can  be  0.974  (assuming  a  density  of 
1026  Kg/m")  . 

The  acceleration  channels  were  calibrated  using  the  fact 
that  when  the  equipment  is  fixed  the  mean  observed  value 
should  be  the  gravitational  acceleration.  The  package  was 
fixed  and  rotated  from  -20°  to  +20°,  in  5°  intervals  on  two 
quasi-orthogonal  planes,  with  one  of  the  planes  containing  one 
of  the  horizontal  accelerometers.  Least  square  fits  to 
straight  lines  were  performed,  resulting  in  the  following 
calibration: 

ACC^=  0 . 00108549  .  .  0324 

ACCy=  0 . 00119582  .  CO£WTSy-19 .967  3 
ACC-=-0 . 00110622  ,  COUNTS^+IS  .9346 


where  ACCi  is  the  acceleration  in  m/s^  and  COUNTSj  the  outputed 
quantities.  The  observed  mean  square  errors  are  0.0532  m/s^, 
0.0210  m/s^  and  0.0634  m/s^  in  the  X,  Y  and  Z  channels 
respectively. 

To  evaluate  the  performance  of  the  equipment,  two  5 
minutes  data  sets  were  selected  from  the  WHARF  #2  Experiment 
at  the  August  2  1991.  This  experiment  will  be  described  in 
more  detail  in  the  next  chapter.  The  calibrated  temperature 
values  are  very  closed  to  those  observed  using  a  mercury 
thermometer  bucket  temperature.  The  accelerometer  mean  values 
are  reasonably  close  to  the  expected.  Because  there  was  no  in 
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situ  calibration,  only  mean  quantities  can  be  analyzed.  A  more 
detailed  description  of  these  data  will  be  presented  in  the 
next  chapter.  This  experiment  was  performed  with  the 
instrument  fixed  to  the  bottom,  so  no  shear  probe  data  was 
analyzed.  No  pressure  data  was  available  from  this  deployment. 

The  first  step  in  the  analysis  is  to  look  at  the  pre- 
processed  doppler  data.  A  sample  of  two  bistatic  reciprocal 
channels  (frequency,  power  and  bandwidth),  and  calibrated 
temperature  data  is  shown  in  the  figures  4.F.1-2. 

A  careful  analysis  was  required  for  the  velocity 
estimation,  because  there  was  no  previous  experience  with  this 
method.  It  was  necessary  to  validate  all  the  processing  steps 
(both  electronic  and  software). 

The  reciprocal  channel  doppler  frequency  estimates  are 
generally  similar,  see  figure  4.F.1,  except  for  individual 
spikes,  or  in  regions  of  very  low  signal  to  noise  ratio  (SNR) 
which  are  also  seen  by  a  large  bandwidth.  Some  of  these 
discrepancies  will  be  attenuated  by  using  a  variable  pulse 
repetition  interval  in  future  deployments,  eliminating 
coherent  additive  interference. 

To  increase  the  confidence  and  deglitch  the  acoustic  data, 
an  adaptative  procedure  has  been  developed.  The  probability 
distribution  function  (PDF)  of  the  differences  histogram  (see 
in  the  Figure  4.F.3),  shows  that  the  differences  between  two 
consecutive  values  are  smaller  than  70  Hz  (0.01  m/s)  75%  of 
the  time.  Also,  the  scatter  plots  between  the  first 
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Fiqurp  4.F.1  -  Reciprocal  Bistatic  channels,  doppler 

frequencies,  received  powers  and  bandwidth. 
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Figure  4.F.2  -  Calibrated  Temperature  data.  Simultaneous  to 
the  data  shown  in  figure  4.F.1  . 

differences  and  the  received  powers  (Figure  4.F.4)  shows  that 
the  bigger  differences  are  mainly  associated  with  low  received 
powers,  corresponding  to  the  cases  when  noise  from  additive 
interference  and  system  noise  may  be  significant.  We  have  two 
independent  realizations  of  the  same  random  variable 
(reciprocal  channels)  that  can  be  used  either  to  increase  the 
confidence  of  the  estimates  or  to  validate  them. 

To  deglitch  the  data  and  obtain  a  clean  averaged  time 
series  of  the  doppler  frequencies,  an  intelligent  adaptative 
algorithm  must  be  used  based  on  the  surrounding  points.  A 
procedure  was  implemented  to  compare  the  signal  at  each  time 
with  a  predicted  value  obtained  from  a  4’^  order  linear 
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differences  histogram. 


predictor  based  on  the  previous  20  point  autocorrelation 
function.  Whenever  the  difference  was  greater  than  four  times 
the  standard  deviation  (also  estimated  from  the  previous  20 
points),  that  point  was  assumed  non-valid.  When  both  channels 
were  assumed  bad,  or  if  a  first  difference  fixed  criteria  was 
violated,  the  predicted  value  was  used.  We  can  expect  an 
improvement  in  the  SNR  of  the  estimates  by  using  this  method. 
Some  examples  of  how  this  procedure  works,  are  shown  in  the 
figure  4 . F . 5 . 

To  transform  these  deglitched  doppler  shifts  to  velocity 
components  in  an  orthogonal  frame  fixed  to  the  instrument,  a 
linear  transformation  must  be  performed.  This  transformation 
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between  first  differences  and  received  power. 

will  not  distribute  the  velocity  variance  evenly  through  the 
3  orthogonal  components  of  velocity,  giving  more  weight  to  the 
component  along  the  instrument  vertical;  this  componeTit  will 
also  have  a  greater  statistical  confidence. 
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DEGLITCHING  FUNCFION  -  u;  ing  mean  or  predicted  values 
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DEGLITCHING  FUNCTION  -  thresholds  and  errors 
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Figure  4.F.5  -  Bl‘3tatic  channels  deglitching  function. 
Whenever  the  difference  curves  go  higher  than  the  threshold 
curves,  in  the  lower  plot,  that  point  was  assumed  as  bad  data. 
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V.  •nJRBULENCE  OBSERVATIONS  USING  THE  CDV  PACKAGE 

A.  WHARF  #2  EXPERIMENT 

The  purpose  of  this  experiment  was  to  asses  the 
performance  of  the  CDV  package  in  oceanic  backscatter 
conditions,  test  its  ability  to  measure  ocean  fluxes  and 
surface  wave  fields,  and  to  study  the  forcing  mechanisms  in 
ocean  surface  layers. 

The  CDV  package  was  deployed  off  the  Monterey  Wharf  #2  on 
2  and  3  August  1991.  The  tripod-mounted  CDV  was  placed  on  the 
sand-bed  with  the  top  of  the  instrument  at  depths  ranging  from 
0  to  0.75  meter.  The  depth  changed  with  the  tide,  and  the 
total  water  depths  ranged  between  1.7  to  2.5  meter. 

The  package  was  controlled  and  monitored  from  a  van  parked 
on  the  wharf,  from  which  observations  of  air  temperature 
(moist  and  dry),  bucket  sea  surface  temperature,  short  and 
long  wave  radiation  and  wind  direction  and  speed.  The 
equipment  deployment  is  sketched  on  figure  5.A.I. 

The  air  temperature  was  measured  using  a  mercury 
psychrometer  and  the  sea  surface  temperature  using  a  protected 
mercury  thermometer.  The  wind  speed  and  direction  was  recorded 
using  a  MRI  anemometer  and  wind  vane,  while  the  short  and  long 
wave  radiation  was  obtained  from  an  Epply  radiometer. 
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Figure  5.A.1  -  Wharf  §2  Experiment  set-up 


The  environmental  conditions  were  characterized  by  an 
overcast  sky,  moderate  to  weak  winds  from  WNW.  Small  amplitude 
swell  from  the  North  was  being  reflected  from  the  solid  walls 
below  the  wharf.  Due  to  the  barrier  created  by  the  wharf  walls 
on  the  wind  field,  the  measured  winds  did  not  correspond  well 
to  those  over  the  CDV  package.  These  light  wind  conditions 
were  not  ideal  to  evaluate  the  overall  performance  of  the 
instrument,  because  we  would  not  expect  to  find  high 
turbulence  levels  in  the  high  frequency  range.  Nevertheless, 
because  we  understand  the  physics  of  the  linear  wave  fields 
and  we  do  not  have  lateral  boundary  effects  (as  we  would 
encounter  in  a  tank  experiment) ,  we  were  able  to  measure  the 
performance  of  the  package  at  the  low  and  medium  frequency 
range  using  the  wave  field  signal,  and  to  estimate  the  noise 
floor  at  higher  frequencies. 

The  CDV  data  were  recorded  on  VCR  tapes  ( labeled  WCDVl  to 
WCDV7),  from  1100  until  2300  on  2  Aug  1991,  covering  the 
afternoon  high  tide  at  1545  pm,  and  from  0200  to  0800  at  3 
Aug,  covering  the  morning  high  tide  at  0515.  The  short  and 
long  wave  radiation  and  the  wind  data  were  recorded 
continuously  from  1100  am  of  2  Aug  until  0800  am  of  3  Aug, 
sampled  at  8  Hz  and  spooled  by  an  HP  3  32  computer  to  an 
optical  drive  for  further  processing.  The  air  and  bucket 
temperature  were  recorded  every  hour  from  1600  2  Aug  to  0800 
3  Aug. 
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After  a  careful  analysis,  only  the  afternoon  CDV  data  from 
1300  to  1800  show  high  data  quality,  while  the  morning  data 
show  large  data  gaps  due  to  the  fact  that  the  equipment  was 
only  partially  submerged. 

Following  DSP  preprocessing,  all  the  data  were  spooled  to 
an  HP  VECTRA  386  and  the  CDV  bistatic  data  to  a  SUN  SPARC 
STATION  for  further  processing. 

After  calibration,  the  8  Hz  sampled  wind  and  radiation 
data  were  10  min  averaged  and  decimated  at  a  1/min  sampling 
rate.  The  CDV  sea  temperature  was  calibrated  and  used  for 
cross-spectral  analysis  and  computation  of  fluxes  over  20 
minutes  averages.  The  bistatic  doppler  frequency  estimates 
were  degl itched  using  the  adaptative  procedure  described  in 
the  previous  chapter,  and  used  to  estimate  the  velocity  vector 
in  an  instrumental  reference  frame.  The  instrumental 
velocities  were  then  properly  rotated  to  a  vertical  reference 
frame.  All  the  processed  data  (intermediate  and  final), 
spectral  and  flux  estimates  and  output  plots  were  stored  on  an 
optical  disk. 

B.  METEOROLOGICAL  DATA 

This  section  describes  and  analyzes  the  tide,  wind, 
radiation,  air  and  sea  surface  temperature  data.  It  will  focus 
only  on  the  period  from  1100  to  2000  of  2  Aug  1991, 
simultaneous  with  the  useable  CDV  data. 
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As  is  shown  in  the  Figure  5.B.1,  within  this  period  there 
was  a  high  tide  at  1545  pm.  The  tide  height  determines  the 
depth  of  the  CDV  measured  quantities;  so  in  Figure  5.B.1,  the 
expected  mean  depth  of  the  top  of  the  CDV  is  shown,  assuming 
that  at  1100  it  was  at  a  0  m  depth.  Also  shown  as  vertical 
straight  lines  are  the  times  of  the  evaluated  ocean  fluxes. 

Due  to  problems  in  the  radiation  sensors,  long  wave 
radiation  data  was  not  available.  This  fact  restricts  our 
ability  to  compute  the  surface  buoyancy  fluxes;  however,  we 
can  use  the  short  wave  radiation  and  wind  stress  as  indicators 


CDV  package  at  2  August  1991 . 
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of  the  forcing  mechanisms  in  the  ocean  upper  layer.  The  wind 
data  does  not  accurately  reflect  wind  stress  at  the  CDV 
location  due  to  the  barrier  created  by  the  wharf  walls,  but  it 
will  show  the  maximum  expected  mechanical  forcing  of  the  ocean 
surface.  The  time  evolution  of  these  variables  is  shown  in 
Figure  5.B.2,  for  the  period  from  1200  to  2000  pm  at  2  Aug 
1991,  based  on  measurements  made  from  the  top  of  the  wharf. 

These  plots  show  that  the  wind  stress  was  decreasing 
slightly  after  1500.  Due  to  the  overcast  conditions  during 
almost  all  the  period  of  the  experiment,  we  should  not  expect 
a  large  variability  on  the  long  wave  radiation,  so  the  changes 
in  the  surface  buoyancy  flux  should  be  determined  primarily  by 
the  short  wave  radiation  and  by  the  surface  layer  stability. 
The  short  wave  radiation  followed  a  cloud  covered  diurnal 
cycle  increasing  until  1500,  and  decreased  thereafter.  The 
peak  shown  at  approximately  1420  corresponds  to  a  short  period 
of  time  of  partially  clear  sky.  From  these  conditions,  we  may 
expect  a  decrease  in  the  weak  mechanical  steering  due  to 
decreasing  winds  and  a  reduction  in  the  stability  of  the  upper 
layer  and  due  to  a  reduction  of  solar  heating.  Therefore,  we 
should  not  expect  to  find  high  levels  of  turbulent  energy  in 
the  form  of  scalar  fluxes  and  Reynolds  stresses,  i.e.  we 
should  not  expect  to  observe  an  active  mixing  layer,  so  that 
the  ocean  conditions  should  be  dominated  by  surface  gravity 
wave  fields. 
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Figure  5.B.2  -  Wharf  #2  Experiment:  Wind  stress  and  short 
wave  radiation,  from  1200  to  2000  at  2  Aug  1991. 


C ,  CDV  DATA 

The  data  between  1500  and  1800  of  2  Aug  1991  were  analyzed 
in  more  detail  and  the  ocean  fluxes  computed  using  an  eddy 
flux  correlation  method  [Ref. 2].  The  first  step  was  to  select 
20  minutes  ensembles  every  1  hour.  These  data  sets  were  then 
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spooled  to  an  optical  drive,  and  the  data  analyzed  in  more 
detail  using  C  programs  and  MATLAB  software. 

After  conversion  to  engineering  units,  the  acceleration 
and  temperature  data  were  averaged  over  all  the  20  minute 
ensembles;  the  tilt  angles  were  computed  and  the  mean 
temperatures  compared  with  the  observed  bucket  temperature. 
TABLE  6  -  CDV  AVERAGED  DATA 


Hour  ( 2  Aug ) 

15:30 

16:30 

17:30 

18:30 

Acceleration-X 

(m/s=) 

0.1626 

0.1588 

0.1630 

0.1553 

Acceleration-y 

(m/s*) 

-0.6645 

-0.6647 

-0.6718 

-0.6787 

Acceleration-z 

(m/s*) 

-9.8278 

-9.8330 

-9.8277 

-9.8288 

Direction-Deg 

(ref. to  x-axis) 

-76 . 25° 

-76.56° 

-76.36° 

-77.10° 

Tilt  (Deg) 

(from  vertical) 

3.98° 

3.976° 

4.02° 

4.05° 

Mean  Temperature 

(bucket) 

17.1°C 

16.8°C 

16.3°C 

16.0°C 

Mean  Temperature 

(CDV  sensor) 

16.24°C 

16.23°C 

16.17°C 

15.82°C 
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Hour  (2  Aug) 

15:30 

16:30 

17:30 

18:30 

CDV  Depth 

0.74m 

0.71m 

0 . 64m 

0 . 36m 

(Tide  Table) 

We  can  see  that  over  the  full  deployment  there  was  a 
constant  tilt  of  4  degrees  from  the  vertical,  and  an  offset  of 
-76.6  degrees  from  the  x-axis,  defined  as  the  bisection  of 
transducers  1  and  2.  These  values  were  used  to  rotate  the 
velocity  data  to  obtain  a  vertical  reference  frame.  Examples 
of  the  measured  velocities  are  shown  in  Figure  5.C.I.  The 
small  magnitude  and  variability  of  these  values  is  consistent 
with  the  observed  conditions  of  low  or  nonexistent  turbulence 
and  small  amplitude  wave  field. 

The  mean  temperature  compares  well  with  the  bucket 
temperature,  but  the  44  H2  measured  temperature,  shows  some 
high  gradients  (as  is  shown  in  Figure  5.C.2)  not  correlated 
with  the  vertical  velocities,  suggesting  high  horizontal 
stratification.  The  strong  and  sharp  variations  of 
temperature,  by  itself,  suggest  horizontal  advection  of  warm 
patches  of  water.  This  hypothesis  is  further  discussed  in  the 
next  section,  when  the  temperature  fluxes  are  discussed. 

D.  CXIEAN  FLUXES  AND  STRATIFICATION 

Figures  5.D.1  to  4,  show  the  spectral  evolution  of 
temperature  and  of  the  3  velocity  components.  The  freguency 
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density  units  of  all  spectral  quantities  are  such  that  the 
integration  over  the  full  frequency  range  gives  the  total 
variance.  These  spectra  are  characterized  by  a  broad  peak, 
more  evident  in  the  velocity  spectra,  at  frequencies  between 
0.2  and  0.6  Hz.  They  represent  short  period,  surface  gravity 
waves,  with  wavelengths  of  the  order  of  6  m  with  e-folding 


at  44  Hz. 

depths  of  the  order  of  1  m  (i.e.  deep  water  waves).  This  fact 
is  consistent  with  the  observed  similar  energy  content  in  the 
vertical  and  horizontal  velocity  components  (10‘*  (m/s)^/Hz). 
The  broad  band  shape  of  the  spectrum  suggests  that  this  energy 
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is  due  to  wind  waves  starting  to  shoal,  probably  generated  in 
the  near  field. 

There  are  also  three  peaks,  possibly  associated  with  low 
frequency  swell  or  seiche,  at  the  frequencies  of:  2.5x10'^  HZ, 
4.5x10'  Hz  and  7x10'^  Hz:  these  spectral  peaks  can  be  seen  in 
the  figure  5.D.5,  which  shows  the  high  resolution  spectra  of 
the  variables  at  1530.  These  peaks  are  more  noticeable  in  the 
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Figure  5.C.2  -  Wharf  ^2  Experiment:  Temperature  data. 


sampled  at  44  Hz. 

spectra  of  the  horizontal  velocity,  significant  to  a  80% 
confidence  level,  and  are  highly  attenuated  in  the  vertical 
velocity  spectra  as  we  can  expect  for  a  shallow  water  wave. 
These  harmonically  related  frequencies  suggest  weak  non-linear 
interactions  (resonance),  with  energy  being  transferred, 
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Figure  5.D.5  -  Wharf  #2  Experiment  -  High  resolution  spectra 
at  15:30, 
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between  the  different  waves  components.  If  ressonance  is 
occuring,  these  components  will  contribute  to  the  presence  of 
stresses  in  the  upper  layer. 

Besides  telling  us  that  the  velocity  field  was  dominated 
by  the  surface  gravity  waves,  the  low  energy  levels  at  high 
frequencies  also  shows  the  noise  floor  of  the  instrument  with 
energy  density  on  the  order  of  10“’  (m/s)^/Hz  in  the  velocity 
spectra,  and  10‘^°  (Deg  C)^/Hz  in  the  temperature  spectra.  This 
value  represents  noise  significant  velocities  of  the  order  of 
0.005  m/s,  that  are  smaller  than  typical  velocities  associated 
with  active  oceanic  turbulence,  typically  of  the  order  of  0.01 
m/s  [Ref. 22].  Therefore,  for  turbulent  flux  observations  we 
should  expect  to  obtain  good  SNR  performances  with  this 
instrument . 

Integrating  the  velocity  spectra,  we  obtain  the  kinetic 
energy  for  each  ensemble.  These  values  and  their  evolution 
over  the  4  analyzed  hours  are  important  in  the 
characterization  of  the  observed  fields.  We  can  see  that  the 
total  horizontal  kinetic  energy  was  significantly  larger  than 
the  vertical,  due  to  the  nature  of  the  wave  fields. 
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HOURS 

15:30 

16:30 

17:30 

18:30 

u-component 

(mVs^) 

0.0044 

0.0049 

0.0062 

0.0119 

v-component 

(mVs") 

0.0107 

0.0118 

0.0149 

0.0238 

In  order  to  compute  and  interpret  the  fluxes  of 
temperature  and  the  Reynolds  stresses,  the  co-spectra,  phase 
spectra  and  coherence  were  computed  for  the  several 
combinations  of  variables.  Using  Parseval  theorem  [Ref. 14], 
integrating  the  co-spectra  over  frequency  provides  estimates 
of  the  scalar  fluxes  or  the  Reynolds  stresses.  Therefore,  the 
co-spectrum  will  tell  us  how  these  fluxes  are  distributed  over 
frequency.  In  figures  5.D.6,  7. a,  7.b  and  7.c,  we  can  see  the 
evolution  of  the  co-spectrum  between  the  perturbed  quantities 
of  the  temperature  (t')  and  vertical  velocity  component  (w') 
hereafter  designated  by  t'w',  and  between  the  three  velocity 
components  u'w' ,  v'w'  and  u'v'.  The  low  values  in  the  co¬ 
spectrum  and  the  negligible  energy  above  1  Hz  are 
characteristic  of  a  linear  wave  field  with  low  turbulence 
levels . 

The  negative  peak  found  at  the  very  low  frequency  range  in 
the  co-spectrum  of  the  horizontal  velocity  components,  Figure 
5.D.7.a,  confirms  the  presence  of  non-linear  interactions 
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between  harmonically  related  wave  components  [Ref. 23],  as 
suggested  by  the  velocity  spectra  in  Figure  5.D.5. 

The  co-spectra  of  temperature  and  vertical  component  of 
the  velocity,  shows  that  the  thermal  energy  flux  is  negative 
in  the  short  period  wave  range  (deep  water  waves),  and 
positive  over  the  long  period  waves  range  at  the  1530  and  1630 
hours.  Later  energy  starts  cascading  into  the  smaller  scales, 
where  previously  the  co-spectral  density  was  almost  zero 
(Figure  5.D.6).  This  is  partially  consistent  with  an  increase 
in  the  horizontal  momentum  flux,  though  it  remains  restricted 
in  the  low  frequency  range.  No  significant  changes  are  seen  in 
the  vertical  momentum  flux.  This  picture  suggests  that  after 
the  1730,  due  to  the  shallow  depth  of  the  instrument  and  to 
the  reduction  of  insolation  generated  stratification, 
turbulent  levels  tend  to  increase. 

After  integrating  equivalent  estimates  for  each  data  set, 
the  time  evolution  of  the  temperature  fluxes  and  stresses  were 
calculated ; 


TABLE  8  -  TEMPERATURE  AND  MOMENTUM  FLUXES 


HOURS 

15:30 

16:30 

17:30 

18:30 

<u'v'> 

(mVs") 

-0.0046 

-0.0051 

-0.0063 

-0.0108 

<u'w'> 

(mVs") 

-0.0001 

-0.0000 

0.0002 

0.0015 
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HOURS 

15:30 

16:30 

17:30 

18:30 

<v' w' > 

(mVs") 

0.0011 

0.0014 

0.0020 

0.0012 

<t 'w'> 

-2.2277 

-1.0943 

-2.5964 

-10.5082 

(ra.°C/s) 

xlO-^ 

xlO'^ 

XlO'^ 

xlO'^ 

This  table  shows  that  the  Reynolds  stresses  and  temperature 
fluxes  slightly  increased  with  time,  as  stratification  was 
decreasing . 

The  phase  spectra  and  coherence  are  tools  that  help  the 
interpretation  of  the  co-spectrum.  The  phase  spectra  quantify 
the  phase  lags  between  variables,  showing  departures  from  the 
in-phase  components  shown  in  the  co-spectra.  These  spectra  are 
useful  in  the  identification  of  quadrature  conditions, 
characterizing  linear  wave  components,  and  in  the  analysis  of 
non-linear  characteristics.  The  coherence  spectrum  shows 
whether  cross-spectral  components  have  constant  phase 
relationships  in  time.  For  linear  wave  fields,  coherence 
should  be  high,  near  1  (no  fluxes  or  stresses),  while  for 
turbulent  or  non-linear  fields  the  coherence  should  have  low 
values.  Examples  of  these  estimates  are  shown  in  the  figures 
5.D.8  and  9  for  the  temperature  vertical  flux  and  Reynolds 
stresses  at  1530.  The  features  shown  in  these  figures  are 
similar  to  the  characteristics  of  the  other  data  sets,  with 
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quadrature  phase  and  high  coherence  through  the  spectral  peaks 
range,  as  expected  for  a  linear  wave  field. 

The  low  correlation  between  temperature  and  vertical 
velocity  hint  that  we  cannot  explain  the  observed  temperature 
variability  only  from  the  vertical  stratification,  in  which 
water  is  pumped  up  and  down  past  the  sensor  following  the  wave 
orbital  motion.  On  other  hand,  if  we  assume  that  horizontal 
stratification  exists,  then  the  horizontal  velocities  may 
explain  this  variability,  with  patches  of  warmer  or  colder 
water  being  advected  by  the  mean  horizontal  flow  and  wave 
orbits.  To  test  this  hypothesis,  the  horizontal  temperature 
fluxes  were  computed  as  is  shown  in  the  figure  5.D.10  for  the 
1530  data  set.  In  fact,  they  do  not  show  a  stronger 
correlation;  however,  the  probable  high  variability  on  the 
horizontal  density  gradient  will  introduce  an  additional 
effect  on  these  estimates.  These  estimates  may  also  be 
misleading  if  not  computed  from  sufficiently  long  averaging 
intervals  as  the  sensor  was  stationary  in  the  field. 

The  coherence  spectra(  Figure  5.D.9)  again  reinforce  the 
fact  that  the  wave  regime  at  the  middle  frequencies  has  a 
linear  character,  showing  high  coherence  between  v  and  w 
within  this  range.  These  characteristics  are  maintained 
throughout  the  observation  period.  Similarly,  the  coherence 
between  the  horizontal  components  of  velocity  is  more 
significant  in  the  low  frequency  wave  range,  but  with  slightly 
lower  values. 
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Figure  5.D.10  -  Wharf  ^2  Experiment  -  Horizontal  fluxes  of 
temperature  at  15:30. 
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VI .  CONCLUSIONS 


The  need  to  explicitly  observe  turbulent  processes  arises 
from  a  wide  class  of  problems  in  physical  oceanography,  as 
described  in  the  review  in  chapter  II.  The  CDV  package  was 
developed  to  simultaneously  measure  temperature,  conductivity, 
pressure,  shear  and  the  3  components  of  velocity  at  sub¬ 
centimeter  scales. 

The  unigue  measurement  technigues  used  in  this  instrument 
require  special  data  processing  and  analysis  methods.  The 
three  component  velocity  estimations  rely  on  a  doppler 
frequency  spectral  estimator,  working  in  real-time,  that  must 
use  a  small  number  of  coherently  sampled  data  points.  In  order 
to  choose  the  meet  suitable  a_gorithm,  and  to  assess  its 
physical  limitations,  Monte-Carlo  simulations  were  performed 
using  several  possible  spectral  estimators.  From  these 
simulations,  the  Pulse  Pair  algorithm,  or  equivalently  the 
MUSIC  order  2  algorithm,  showed  the  best  performance,  followed 
by  the  Peak  FFT  algorithm.  The  FFT  estimator,  corresponds  to 
a  maximum  likelihood  estimator,  but  suffers  from  the  fact  that 
for  a  small  number  of  points,  the  discretization  errors  are 
very  dominant,  and  it's  performance  is  seriously  degraded. 
Using  these  results,  the  Pulse  Pair  estimator  was  implemented, 
together  with  the  pre-processing  of  the  other  measured 
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quantities,  on  a  Digital  Signal  Processor  for  real-time  pre¬ 
processing  and  review. 

To  test  the  CDV  software  and  hardware,  the  equipment  was 
deployed  for  the  first  time  off  Monterey  Wharf  #2,  on  the  2 
and  3  of  August  1991,  along  with  measured  environmental 
conditions.  This  first  analysis  showed  that  the  performance  of 
the  bistatic  velocity  estimation  is  below  that  expected,  with 
some  conspicuous  spikes  superimposed  rn  valid  data.  By 
introducing  a  change  in  the  transmitting  cycles  of  the 
transducers  this  problem  should  be  overcome  in  future 
deployments.  The  calibration  of  the  other  scalar  variables  and 
acceleration  gave  good  results.  To  deglitch  the  bistatic 
doppler  frequencies  and  increase  its  confidence,  we  used  the 
redundancy  in  the  velocity  components  observations,  with  two 
independent  estimates  from  each  of  the  two  reciprocal  acoustic 
bistatic  paths.  An  adaptative  procedure  was  developed,  relying 
on  a  one  lag  linear-predictor  and  predictively  choosing  the 
best  estimate  of  the  doppler  frequency.  This  procedure 
significantly  improved  the  SNR  of  the  velocity  estimates. 

The  CDV  Wharf  data  shows  a  surface  gravity  wave  field  with 
energy  in  two  different  ranges,  superimposing  near  field 
generated  high  frequency  wind  waves  with  a  low  frequency 
swell.  The  wind  waves  sho>?  the  expected  characteristics  of 
deep  water  waves  while  the  low  frequency  waves  clearly  have 
shallow  water  waves  characteristics.  The  high  resolution 
spectra  and  the  co-spectra  show  that  resonance  interactions 


may  be  occurring  at  the  low  frequency  range.  The  high 
frequency  wind  waves  show  a  stronger  coherence  and  do  not 
contribute  significantly  to  the  momentum  flux,  but 
nevertheless  the  major  contribution  to  the  temperature  flux 
occur  over  these  frequency  ranges. 

The  fact  that  the  observed  temperature  variability  did  not 
correlate  well  with  the  measured  vertical  velocity,  even  at 
the  deep  water  wave  frequency  range,  suggest  that  an  irregular 
horizontal  stratification  existed,  with  pools  of  warm  or  cold 
water  being  advected  past  the  sensor,  by  the  wave  orbits  and 
mean  flow. 

From  the  vertical  velocity  and  temperature  co-spectra 
analysis,  we  also  observed  an  increase  in  the  temperature  flux 
associated  with  a  shallow  depth  and  reduction  in  the  vertical 
stratification,  following  the  decrease  in  the  solar  heating 
after  1730.  However,  due  to  remaining  significant 
stratification  no  significant  turbulence  develops,  as  we  can 
observe  in  the  co-spectra  of  the  velocity  components. 

At  the  high  frequency  ranges,  we  did  not  found  significant 
energy  due  to  the  absence  of  strong  turbulence,  allowing  the 
instrumentation  noise  floor  to  be  seen.  This  noise  floor 
represents  velocities  in  the  order  of  0.005  m/s,  values  below 
the  ones  that  we  should  expect  to  find  during  strong  turbulent 
conditions.  Because  this  data  was  the  first  that  was  analyzed 
with  this  equipment,  we  expect  in  future  deployments  to 
improve  the  SNR  in  the  velocity  measurements  by  an  order  of 
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magnitude  by  reducing  the  acoustic  noise  component  5  to  10 
times  and  altering  acoustic  sampling  methods.  However,  the 
existing  noise  floor  is  already  below  the  minimum  requirements 
for  studying  energetic  turbulence. 

The  CDV  data  seems  consistent  with  the  observed 
environmental  conditions,  confirming  the  expected  capabilities 
of  the  instrument  to  observe  real  data,  including  wind  wave 
fields,  and  showing  the  noise  limitations  that  we  should 
expect  when  observing  turbulence  fields.  However,  in  order  to 
further  study  its  precision,  this  data  should  also  be  compared 
with  simultaneous  observations  of  the  ocean  variables  measured 
by  independent  instruments.  The  CDV  package  will  be  able  to 
perform  simultaneous  measurements  of  oceanographic  variables 
at  the  sub-centimeter  scales,  determining  all  terms  in  the 
Reynolds  equations.  This  should  improve  our  characterization 
of  turbulent  fields,  with  applications  in  a  variety  of 
oceanographic  turbulence  problems.  These  applications  will  be 
explored  using  the  instrument  in  a  loose-tethered  profiling 
mode . 
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